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ABSTRACT 


This  study  examines  various  computational  techniques  to  analyze  dynamic  response  and 
failure  of  sandwich  composite  materials  subject  to  fluid- structure  interaction  character¬ 
ized  by  an  acoustic  field  or  the  propagation  of  velocity  potential  according  to  the  wave 
equation.  A  displacement-only  plate  finite  element  is  developed  and  implemented  using 
Discontinuous  Galerkin  (DG)  methodology;  its  accuracy  compares  favorably  to  both  the¬ 
ory  and  Continuous  Galerkin  methods.  Several  approaches  to  analyzing  debonding  failure 
between  skin  and  core  layers  of  sandwich  composite  structures  are  demonstrated  and  eval¬ 
uated;  partial  disconnection  between  neighboring  elements  at  a  debonding  site  shows  good 
qualitative  agreement  with  known  physical  phenomena.  A  hybrid  Finite  Element-Cellular 
Automata  (FE+CA)  approach  to  modeling  an  acoustic  field  with  non-reflecting  boundary 
conditions  is  presented,  validated  numerically,  and  favorably  compared  with  experimen¬ 
tal  results.  The  FE+CA  fluid  model  is  then  combined  with  the  DG  structural  model  to 
simulate  fluid- structure  interaction;  this  combined  model  compared  favorably  with  experi¬ 
mental  results  for  the  strain  field  of  laminated  plates  subject  to  low-velocity  impact.  Each 
technique  addressed  shows  promise  for  flexible  and  accurate  modeling  of  failure  initiation 
and  propagation  in  sandwich  and  laminate  composites  subject  to  fluid-structure  interaction 
with  moderate  computational  costs. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

A.  MOTIVATION  AND  OBJECTIVE .  1 

B.  PRIOR  WORK  .  2 

C.  ORGANIZATION .  4 

II.  STRUCTURAL  MODEL .  7 

A.  LINEAR  ELASTICITY .  7 

B.  DISCONTINUOUS  GALERKIN  FORMULATION .  9 

1.  Solid  Element  .  9 

a.  Validation .  19 

2.  Plate  Element  .  19 

a.  Validation .  21 

3.  Stabilization  and  Penalty  Parameters .  25 

III.  COMPOSITE  STRUCTURES .  27 

A.  MULTI-SCALE  MODEL .  27 

B.  SANDWICH  COMPOSITES .  28 

C.  FAILURE  MODES  AND  IDENTIFICATION .  32 

D.  FAILURE  MODELING .  37 

1.  Damage  via  Complete  Disconnection .  37 

2.  Damage  via  Partial  Disconnection .  42 

3.  Damage  via  Reduced  Moduli .  47 

E.  SYNTHESIS  50 

IV.  FLUID  MODEL .  51 

A.  VELOCITY  POTENTIAL  AND  WAVE  EQUATION  THEORY  .  51 

B.  FINITE  ELEMENT  MODEL  OF  THE  WAVE  EQUATION  ....  53 

C.  CELLULAR  AUTOMATA  MODEL  OF  THE  WAVE  EQUATION  54 

1.  One  Dimension .  55 

2.  Boundary  Conditions .  57 

3.  Discretization  and  Model  Fidelity  .  58 

4.  Convergence .  62 

5.  Three  Dimensions .  63 

vii 


D.  COUPLING  OF  FINITE  ELEMENT  AND  CELLULAR  AUTOMATA 

MODELS  OF  THE  WAVE  EQUATION .  68 

1.  Comparison  with  Homogeneous  Fluid  Domain . 73 

V.  RESULTS .  75 

A.  ACOUSTIC  FIELD  FLUID-STRUCTURE  INTERACTION  ...  75 

1.  Information  Exchange .  75 

2.  Homogeneous  Isotropic  Single-Layer  Structure . 75 

3.  Two-layer  Plates .  78 

4.  Three-layer  Plates .  81 

B.  EXPERIMENTAL  VALIDATION .  86 

VI.  CONCLUSION  .  95 

A.  SUMMARY  OF  FINDINGS  .  95 

B.  FUTURE  WORK .  96 

APPENDIX  A.  TIME  INTEGRATION  ALGORITHMS .  98 

A.  NEWMARK-/3  METHOD  .  99 

B.  a-METHOD . 100 

C.  TIME  DISCONTINUOUS  GALERKIN  METHOD . 100 

APPENDIX  B.  IMPLEMENTATION  DETAILS . 103 

A.  MESHING . 103 

B.  NUMERICAL  INTEGRATION . 103 

C.  APPLICATION  OF  BOUNDARY  CONDITIONS  AND  EXTER¬ 
NAL  LOADS . 104 

D.  MATLAB  SPECIFICS  . 105 

1.  Sparse  Matrices . 105 

2.  CA  Implementation . 105 

E.  MATERIAL  PROPERTIES . 106 

LIST  OF  REFERENCES . 109 


INITIAL  DISTRIBUTION  LIST 


115 


LIST  OF  FIGURES 


Figure  1 .  Canonical  hexahedral  finite  element .  9 

Figure  2.  Connectivity  between  adjacent  Continuous  Galerkin  (CG)  elements. 

After  [31] .  16 

Figure  3.  Connectivity  between  adjacent  Discontinuous  Galerkin  (DG)  elements. 

After  [31] .  17 

Figure  4.  Sparsity  pattern  of  CG  global  stiffness  matrix  for  a  4x2x1  element 

structure .  17 

Figure  5.  Sparsity  pattern  of  DG  global  stiffness  matrix  for  a  4x2x1  element 

structure .  18 

Figure  6.  Convergence  of  three-dimensional  DG  formulation  vs.  discretization 

for  cantilever  beam  deflection .  19 

Figure  7.  Convergence  of  three-dimensional  DG  formulation  vs.  discretization 

for  simply-supported  beam  deflection  .  20 

Figure  8.  Convergence  of  center  deflection  of  a  clamped  plate  comprised  of 

plate  elements .  22 

Figure  9.  Comparison  of  displacement  calculated  for  Continuous  Galerkin  and 

Discontinuous  Galerkin  models  of  a  center  loaded  clamped  plate  ...  23 
Figure  10.  Comparison  of  velocity  calculated  for  Continuous  Galerkin  and  Dis¬ 
continuous  Galerkin  models  of  a  center  loaded  clamped  plate .  23 

Figure  11.  Convergence  of  central  deflection  of  uniformly  loaded  clamped  lam¬ 
inated  plate,  comparison  with  [34] .  24 

Figure  12.  Multi-scale  analysis  cycle  for  a  fibrous  composite.  From  [40] .  28 

Figure  13.  Three  and  Five  Layer  Sandwich  Plate  models .  29 

Figure  14.  Convergence  of  DG  Sandwich  Plate  max  deflection  (due  to  uniform 

load)  to  predicted,  comparison  with  [41,  42] .  31 

Figure  15.  Relative  convergence  of  DG  Sandwich  Plate  max  deflection  and  skin 

stress  (due  to  concentrated  load) .  31 

Figure  16.  Planar  stress  on  the  skin  of  a  clamped  three-layer  sandwich  plate  sub¬ 
ject  to  concentrated  center  force  .  33 

Figure  17.  Planar  stress  on  the  core  of  a  clamped  three-layer  sandwich  plate  sub¬ 
ject  to  concentrated  center  force  .  34 

Figure  18.  Planar  stress  on  the  skin  of  a  clamped  five-layer  sandwich  plate  sub¬ 
ject  to  concentrated  center  force  .  34 

Figure  19.  Planar  stress  on  the  core  of  a  clamped  five-layer  sandwich  plate  sub¬ 
ject  to  concentrated  center  force  .  35 


IX 


Figure  20.  Planar  stress  on  the  bottom  of  a  resin  layer  of  a  clamped  five-layer 

sandwich  plate  subject  to  concentrated  center  force .  35 

Figure  21.  Planar  stress  on  the  top  of  a  resin  layer  of  a  clamped  five-layer  sand¬ 
wich  plate  subject  to  concentrated  center  force .  36 

Figure  22.  Planar  stress  on  the  skin  of  a  damaged  via  complete  disconnection 

clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force  39 

Figure  23.  Planar  stress  on  the  core  of  a  damaged  (via  complete  disconnection) 

clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force  39 

Figure  24.  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (via  com¬ 
plete  disconnection)  clamped  five-layer  sandwich  plate  subject  to  con¬ 
centrated  center  force .  40 

Figure  25.  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (via  complete 
disconnection)  clamped  five-layer  sandwich  plate  subject  to  concen¬ 
trated  center  force .  40 

Figure  26.  Static  deflection  along  centerline  of  damaged  (via  complete  discon¬ 
nection)  clamped  five-layer  sandwich  plate  subject  to  concentrated 
center  force .  41 

Figure  27.  Planar  stress  on  the  skin  of  a  damaged  (by  partial  disconnection  clamped) 

five-layer  sandwich  plate  subject  to  concentrated  center  force . 43 

Figure  28.  Planar  stress  on  the  core  of  a  damaged  (by  partial  disconnection  clamped) 

five-layer  sandwich  plate  subject  to  concentrated  center  force . 43 

Figure  29.  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (by  partial 
disconnection)  clamped  five-layer  sandwich  plate  subject  to  concen¬ 
trated  center  force .  44 

Figure  30.  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (by  partial 
disconnection)  clamped  five-layer  sandwich  plate  subject  to  concen¬ 
trated  center  force .  44 

Figure  31.  Static  deflection  along  centerline  of  damaged  (via  partial  disconnec¬ 
tion)  five-layer  sandwich  plate  subject  to  concentrated  center  force  .  .  45 
Figure  32.  Planar  stress  on  the  skin  of  a  damaged  (by  partial  disconnection) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force, 

fine  view .  45 

Figure  33.  Planar  stress  on  the  core  of  a  damaged  (by  partial  disconnection) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force, 

fine  view .  46 

Figure  34.  Planar  stress  in  the  resin  layer  of  a  damaged  (by  partial  disconnec¬ 
tion)  clamped  five-layer  sandwich  plate  subject  to  concentrated  center 
force,  fine  view .  46 


x 


Figure  35.  Planar  stress  on  the  skin  of  a  damaged  (by  reduced  modulus)  clamped 

five-layer  sandwich  plate  subject  to  concentrated  center  force . 47 

Figure  36.  Planar  stress  on  the  core  of  a  damaged  (by  reduced  modulus)  clamped 

five-layer  sandwich  plate  subject  to  concentrated  center  force . 48 

Figure  37.  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (by  re¬ 
duced  modulus)  clamped  five-layer  sandwich  plate  subject  to  con¬ 
centrated  center  force .  48 

Figure  38.  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (by  reduced 
modulus)  clamped  five-layer  sandwich  plate  subject  to  concentrated 

center  force .  49 

Figure  39.  Center  node  (black)  and  its  Von  Neumann  neighbors  (white) .  55 

Figure  40.  Initial  perturbation  of  an  infinite  string .  56 

Figure  41.  CA  solution  vs.  D’Alembert’s  solution  to  the  one-dimensional  wave 

equation  in  an  infinite  string  .  57 

Figure  42.  Virtual  cell  values  for  various  boundary  conditions  in  one  dimension. 

After  [25] .  58 

Figure  43.  Application  of  various  boundary  conditions  to  CA  calculation  of  one¬ 
dimensional  wave  propagation .  59 

Figure  44.  Initial  perturbations  of  varying  widths .  60 

Figure  45.  CA  solution  to  Id  wave  equation  convergence  as  a  function  of  dx  for 

various  initial  conditions .  60 

Figure  46.  Critical  discretization  vs.  initial  perturbation  width  .  61 

Figure  47.  Convergence  of  CA  wave  equation  rule  to  analytic  solution  as  a  func¬ 
tion  of  da: .  62 

Figure  48.  Three-dimensional  wave  model  at  domain  origin:  time  vs.  iterations  .  64 

Figure  49.  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturba¬ 
tion:  time  vs.  iterations .  65 

Figure  50.  Three-dimensional  wave  model  at  a  point  outside  the  initial  perturba¬ 
tion:  time  vs.  iterations .  65 

Figure  51.  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturba¬ 
tion:  analytic  solution  vs  time,  CA  solution  vs  —  .  66 

Figure  52.  Three-dimensional  wave  model  at  a  domain  origin .  66 

Figure  53.  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturbation  67 

Figure  54.  Three-dimensional  wave  model  at  a  point  outside  the  initial  perturbation  67 

Figure  55.  Finite  Element  fluid  domain  surrounded  by  Cellular  Automata  fluid 

domain  .  69 

Figure  56.  Node  sets  for  exchange  of  data  between  fluid  domains . 70 


xi 


Figure  57.  Comparison  of  velocity  potential  propagation  between  finite  element 
and  cellular  automata  models  with  common  Dirichlet  boundary  con¬ 
ditions  .  71 

Figure  58.  Velocity  potential  propagation  between  finite  element  and  cellular  au¬ 
tomata  domains,  velocity  potential  (0(top))  specified  . 71 

Figure  59.  Velocity  potential  propagation  between  finite  element  and  cellular  au¬ 
tomata  domains,  velocity  (y(top))  specified  .  72 

Figure  60.  Velocity  potential  propagation  between  finite  element  and  cellular 
automata  domains,  FE  inside  CA,  velocity  (v(top))  specified,  non¬ 
reflecting  boundary  conditions .  72 

Figure  61.  Comparison  of  velocity  potential  at  mid-domain  resulting  from  spec¬ 
ified  value  on  one  face:  FE  +  CA  domain  vs.  homogeneous  CA  domain  73 

Figure  62.  Displacement  of  clamped  plate  with  and  without  fluid- structure  inter¬ 
action  .  76 

Figure  63.  Displacement  of  clamped  plate  of  double  modulus  with  and  without 

fluid-structure  interaction .  77 

Figure  64.  Displacement  of  clamped  plate  of  double  density  with  and  without 

fluid-structure  interaction .  77 

Figure  65.  Displacement  of  damaged  clamped  two  layer  E-glass  plate  with  and 

without  fluid-structure  interaction  .  79 

Figure  66.  Velocity  of  damaged  clamped  plate  two  layer  E-glass  plate  with  and 

without  fluid-structure  interaction  .  79 

Figure  67.  Strain  at  center  of  clamped  two  layer  E-glass  plate  with  and  without 

fluid-structure  interaction .  80 

Figure  68.  Strain  at  center  and  edge  of  damage  zone  of  clamped  two  layer  E- 

glass  plate  with  and  without  fluid- structure  interaction .  80 

Figure  69.  Displacement  of  clamped  three  layer  E-glass  plate  with  and  without 

fluid-structure  interaction  and  with  and  without  damage .  81 

Figure  70.  Velocity  of  clamped  three  layer  E-glass  plate  with  and  without  fluid- 

structure  interaction  and  with  and  without  damage .  82 

Figure  71.  Strain  of  clamped  three  layer  E-glass  plate  with  and  without  fluid- 
structure  interaction  and  with  and  without  damage  at  center  of  lower 

E-glass  layer .  82 

Figure  72.  Strain  of  clamped  three  layer  E-glass  plate  with  and  without  fluid- 
structure  interaction  and  with  and  without  damage  at  center  of  inter¬ 
face  layer .  83 

Figure  73.  Strains  in  dry  clamped  three  layer  E-glass  plate  with  and  without 

damage .  84 

xii 


Figure  74.  Strains  in  wet  clamped  three  layer  E-glass  plate  with  and  without 

damage .  85 

Figure  75.  Schematic  of  Vacuum  Assisted  Resin  Transfer  Molding  technique 

(VARTM)  for  plate  manufacture.  From  [49] .  86 

Figure  76.  Drop  Weight  Rig  used  in  Impact  Testing.  From  [48] .  87 

Figure  77.  Drop  Weight  Rig  as  used  for  Impact  Testing  with  Fluid- Structure 

Interaction  (FSI).  From  [48]  88 

Figure  78.  Raw  and  smoothed  experimental  force  data  for  dry  plate  .  89 

Figure  79.  Measured  versus  calculated  strain,  dry  plate,  gage  1 .  90 

Figure  80.  Measured  versus  calculated  strain,  dry  plate,  gage  2 .  90 

Figure  81.  Measured  versus  calculated  strain,  dry  plate,  gage  3 .  91 

Figure  82.  Raw  and  smoothed  experimental  force  data  for  wet  plate .  91 

Figure  83.  Measured  versus  calculated  strain,  wet  plate,  gage  1 .  92 

Figure  84.  Measured  versus  calculated  strain,  wet  plate,  gage  2 .  92 

Figure  85.  Measured  versus  calculated  strain,  wet  plate,  gage  3 .  93 

Figure  86.  Temporal  elements  for  TDG  method.  From  [55] . 102 


xiii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


LIST  OF  ACRONYMS  AND  NOMENCLATURE 


CA  Cellular  Automata 

CG  Continuous  Galerkin 

5P  penalty  parameter 

DG  Discontinuous  Galerkin 

dof  degree(s)  of  freedom 

ex  x  component  of  strain 

£y  y  component  of  strain 

FE  Finite  Element 

FEM  Finite  Element  Method 

FSI  Fluid-Structure  Interaction 

I  IPG  Incomplete  Interior  Penalty  Galerkin 

LDG  Local  Discontinuous  Galerkin 

NIPG  Non-symmetric  Interior  Penalty  Galerkin 

OBB  Oden,  Babuska,  and  Baumann 

SIPG  Symmetric  Interior  Penalty  Galerkin 

ax  x  component  of  stress 

9dg  DG  method  selector 

TDG  Time  Discontinuous  Galerkin 

TSF  time  scaling  factor 

VARTM  Vacuum  Assisted  Resin  Transfer  Molding  technique 


xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xvi 


LIST  OF  TABLES 


Table  1.  Comparison  of  maximum  planar  stress  values  (in  MPa)  and  locations 
for  damaged  and  undamaged  clamped  sandwich  plates  subject  to  a 
concentrated  center  force .  50 

Table  2.  Material  Properties  of  aluminum  honeycomb.  From  [41],  [42]  ....  106 

Table  3.  Material  Properties  of  aluminum  skins,  From  [41],  [42] . 106 

Table  4.  Material  Properties  of  E-glass,  From  [56]  106 

Table  5.  Material  Properties  of  Epoxy  Resin,  From  [57]  107 


xvii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xviii 


ACKNOWLEDGEMENTS 


This  work  is  dedicated  to  the  memory  of  my  grandparents,  Thomas  P.  and  Doris  E.  Lawn 
and  Edward  A.  and  Margaret  F.  Craugh,  who  taught  in  word  and  deed  the  value  of  educa¬ 
tion.  It  would  not  have  been  possible  without  the  love  and  support  of  my  parents,  sisters, 
and  their  families;  in  particular,  my  nieces  Sasha,  Sarah,  and  Anna  Maslowski  and  my 
goddaughter  Sarah  Haynes  inspire  me  constantly. 

I  am  also  indebted  to  Professor  Young  Kwon  for  his  guidance  and  apparently  inex¬ 
haustible  patience. 


xix 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xx 


I.  INTRODUCTION 


A.  MOTIVATION  AND  OBJECTIVE 

Composite  materials  can  have  very  beneficial  properties  for  structural  applications. 
In  particular,  polymer  composites  generally  have  low  density  but  high  strength  and  stiff¬ 
ness,  or  a  very  high  specific  strength  and  stiffness.  Their  resistance  to  corrosion  is  another 
tremendous  advantage.  Consequently,  composite  materials  are  used  in  a  number  of  both 
civil  and  military  applications.  Aerospace  structures  are  among  the  major  applications  of 
polymer  composite  materials,  especially  carbon  fiber  composites.  Increasingly,  composite 
materials  are  used  in  marine  structures,  making  their  use  in  Naval  applications  likewise 
more  common. 

One  of  the  major  difficulties  associated  with  design  and  analysis  of  composite  struc¬ 
tures  is  their  anisotropic  material  properties  and  complex  failure  modes  and  mechanisms 
when  compared  to  metallic  structures.  Anisotropic  composite  material  behavior  can  be  tai¬ 
lored  to  optimize  their  use  in  structures;  however,  diverse  and  inter-connected  failure  modes 
in  composites  remain  a  challenge  for  the  composite  community.  Traditional  metallic  struc¬ 
tures  differ  from  composites  in  one  other  key  aspect:  resilience  and  ductility.  That  is,  their 
ability  to  deform  and  recover  in  the  elastic  zone  gives  the  designer  a  safety  margin  much 
larger  than  that  found  in  the  composite  world.  Composite  materials  have  long  been  used 
as  primary  hull  materials  for  sailboats  and  other  small  craft,  but  the  nature  of  larger  Naval 
applications  of  composite  materials  demands  a  fuller  understanding  of  the  survivability  and 
mission  impacts  of  using  such  structures  in  the  marine  environment. 

A  great  deal  of  research  in  FSI  has  been  conducted  for  aerospace  applications,  tend¬ 
ing  to  focus  on  the  effect  of  the  structure  on  the  fluid  field.  For  marine  applications  with 
a  polymer  composite  structure  in  contact  with  water,  the  comparable  density  of  the  com¬ 
posite  materials  to  that  of  water  results  in  a  significant  hydrodynamic  mass  effect  on  the 
structure.  Therefore,  the  goal  of  this  study  is  to  develop  computational  techniques  to  model 
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and  simulate  transient  dynamic  responses  and  failures  of  composite  structures  with  FSI.  To 
this  end,  it  is  necessary  to  develop  computational  techniques  for  composite  structures  as 
well  as  the  fluid  medium.  Eventually,  a  method  to  couple  both  structural  and  fluid  solvers 
is  developed. 

B.  PRIOR  WORK 

The  Finite  Element  Method  (FEM)  is  a  well  established  tool  for  generating  numer¬ 
ical  solutions  to  the  partial  differential  equations  that  describe  a  wide  range  of  physical 
phenomena.  In  structural  solvers,  the  method  combines  the  geometry  and  constraints  of 
the  structure  with  the  material  properties  of  its  components  to  generate  a  response  (e.g., 
displacement,  stress,  and/or  strain)  to  given  loading.  This  is  accomplished  by  treating  the 
structure  as  a  collection  of  smaller  domains  (finite  elements)  of  relatively  simpler  geome¬ 
try,  applying  and  solving  the  simpler  problem  on  a  smaller  scale  with  constraints  dictated 
by  internal  compatibility,  and  then  combining  the  smaller  solutions  into  a  global  solution. 
This  collection  of  potentially  tedious  but  simple  calculations  is  an  ideal  job  for  a  com¬ 
puter.  The  method  of  weighted  residual  used  in  most  finite  element  codes  is  a  Galerkin 
method  in  which  the  test  function  is  the  derivative  of  the  trial  function  with  respect  to  the 
unknown  variable(s).  Requiring  continuity  between  neighboring  elements  makes  physical 
sense  and  is  the  norm,  but  Discontinuous  Galerkin  (DG)  methods  have  become  an  active 
area  of  research  in  the  last  several  decades. 

DG  approaches  to  solving  boundary-value  problems  have  their  genesis  in  the  work 
of  Nitsche  [1]  who  first  proposed  weak  enforcement  of  boundary  conditions.  Douglas 
and  Dupont  [2],  Arnold  [3],  Baker  [4],  and  Wheeler  [5]  expanded  the  concept  to  weak 
enforcement  of  continuity  between  elements  and  applied  it  to  elliptic  problems.  Easing 
of  the  continuity  constraints  between  elements  suggests  a  number  of  potential  advantages 
to  this  methodology:  element-wise  computations  lend  themselves  to  parallel  computing, 
the  independence  of  the  elements  from  each  other  allows  different  orders  of  interpolation 
within  neighboring  elements,  and  the  potential  to  model  failure  without  having  to  re-mesh 
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is  intriguing.  In  recent  years  discontinuous  approaches  have  been  increasingly  applied  to 
elliptic  partial  differential  equations  such  as  those  that  dominate  linear  elasticity.  Arnold 
and  colleagues  have  produced  two  excellent  overviews  of  the  development  and  characteris¬ 
tics  of  various  methods  in  [6]  and  [7].  In  particular,  they  gathered  the  various  formulations 
and  cast  them  all  in  both  flux  and  primal  forms  to  more  clearly  see  the  differences  and 
identifying  characteristics.  Castillo  [8]  conducted  a  cost  and  performance  analysis  of  three 
of  these  methods.  Brezzi,  Cockbum,  Marini,  and  Suli  [9]  discuss  the  stabilization  mecha¬ 
nisms  necessary  for  effective  DG  formulations.  Specific  work  in  DG  for  solid  mechanics 
can  be  found  in  the  List  of  References,  [10]  -  [22]  among  many  others. 

Not  all  parts  of  a  finite  element  analysis  are  of  equal  interest.  A  homogeneous  beam 
or  plate,  for  example,  can  be  modeled  with  relative  fidelity  with  just  a  few  elements;  while 
the  analysis  of  air-flow  over  a  golf  ball  would  require  elements  of  a  size  comparable  to  the 
ball’s  dimples,  which  could  lead  to  an  exorbitant  computational  cost  for  a  domain  on  the 
order  of  three  ball  diameters.  The  obvious  solution  to  such  a  dilemma  is  to  selectively  refine 
the  mesh  in  the  area  of  interest  while  using  a  mesh  as  coarse  as  possible  in  other  parts  of  the 
domain.  This  work  seeks  to  explore  the  potential  of  using  Discontinuous  Galerkin  (DG) 
finite  element  methodology  to  enable  failure  to  initiate  in  its  natural  location  anywhere  in 
the  domain. 

In  analysis  of  structures  consisting  of  composite  materials,  multiple  scales  of  anal¬ 
ysis  can  be  used  in  conjunction  with  each  other.  That  is,  the  analysis  can  proceed  from 
the  micro-  (atomic  or  molecular)  level  to  the  meso-  (material)  level  and  on  to  the  macro- 
structural)  level  and  back  down  again  as  necessary  with  smearing  or  decomposition  of 
properties  and  loads  as  dictated  by  the  analysis  required.  The  properties  of  a  laminated 
fibrous  composite,  for  example,  are  a  function  of  the  properties  of  both  fiber  and  matrix  as 
well  as  the  weave  pattern  of  individual  laminae  and  lay-up  pattern  of  the  assembled  lami¬ 
nate.  Similarly,  failure  of  the  same  laminate  can  result  from  separation  between  laminae, 
separation  of  fibers  from  matrix,  or  failure  of  individual  fibers. 
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To  accurately  model  FSI  an  appropriate  fluid  model  is  also  required.  In  this  work  no 
attempt  is  made  to  model  or  solve  the  full  Navier-Stokes  equations,  a  subset  that  will  allow 
compatibility  between  fluid  and  solid  regimes  and  approximate  the  hydrodynamic  pressure 
or  acoustic  field  will  suffice.  Olson  and  Bathe  [23]  developed  a  directly  coupled  formula¬ 
tion  that  solves  for  the  velocity  potential  and  hydrostatic  pressure  in  the  fluid  domain  and 
displacements  in  the  structural;  this  will  be  the  starting  point  for  development  of  the  fluid 
model  in  this  work.  In  order  to  apply  the  fluid  model  to  a  maritime  domain,  appropriate 
boundary  conditions  must  also  be  included.  Two  general  approaches  are  to  model  a  vast 
domain  and  concern  oneself  with  a  small  subset  relatively  far  from  simple  but  inaccurate 
boundaries-a  computationally  expensive  proposition-or  to  model  an  appropriately  sized 
domain  with  non-reflecting  boundary  conditions-a  challenging  proposition  that  remains  an 
active  area  of  research  [24].  Application  of  Cellular  Automata  (CA)  to  modeling  the  acous¬ 
tic  field  following  from  the  work  of  Chopard  [25,  26],  Krutar  et  al.  [27],  and  Kwon  and 
Hosoglu  [28]  will  be  explored. 

Increasing  application  of  composite  technology  in  maritime  applications  requires 
further  work  examining  the  mass  effects  imparted  to  composite  structures  in  contact  with  a 
fluid  environment.  In  particular,  evaluation  of  damage  and  residual  strength  are  necessary 
for  a  proper  evaluation  of  the  survivability  potential  of  a  proposed  design. 

C.  ORGANIZATION 

The  balance  of  this  thesis  is  organized  as  follows:  Chapter  II  includes  a  short  re¬ 
view  of  linear  elasticity  followed  by  development  of  both  full  three-dimensional  and  plate 
elements  using  DG  techniques.  Chapter  III  includes  further  discussion  of  the  advantages 
and  challenges  of  composite  materials  and  sandwich  plates,  validation  of  the  developed 
DG  structural  model’s  applicability  to  sandwich  constructs,  followed  with  an  examination 
of  various  failure  models.  Development  and  validation  of  the  fluid  model  incorporating 
velocity  potential  formulation  and  non-reflecting  boundary  conditions  comprises  Chapter 
IV.  Chapter  V  contains  a  demonstration  of  the  coupling  of  the  fluid  and  structural  models 
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as  well  as  comparisons  with  experimental  work  [29].  Conclusions  and  recommendations 
for  continuing  research  in  this  vein  will  be  addressed  in  Chapter  VI. 
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II.  STRUCTURAL  MODEL 


This  chapter  will  discuss  the  formulation  and  various  components  of  the  structural 
models  used  in  this  work.  After  first  reviewing  the  equation(s)  to  be  solved  and  traditional 
numerical  approaches,  Discontinuous  Galerkin  (DG)  methods  will  be  discussed  in  general 
followed  by  a  detailed  formulation  for  a  three-dimensional  solid  element.  That  approach 
will  be  modified  to  model  a  plate  element.  A  discussion  of  the  sensitivities  of  each  element 
to  discretization  and  penalty  parameters  will  ensue. 

A.  LINEAR  ELASTICITY 

In  Voigt  notation,  the  governing  equation  for  the  static  structural  model  is  the  equa¬ 
tion  of  equilibrium, 

V  ■<?+/  =  ()  (1) 

in  which 

®  x  &y  & z  Try  7~yZ  TXz} 

is  the  stress  vector  and 

/=  {/.  fy  A}T  (3) 

is  the  body  force  vector.  Combining  the  constitutive  equation  of  generalized  Hooke’s  Law, 

<?  =  [D]e  (4) 

where  e  =  {sx  ey  ez  ^xy  ^yz  ^XZ}T  is  the  strain  vector  and  \D]  is  the  6x6  material 
property  matrix,  and  the  kinematic  strain-displacement  relationship, 

£  =  ^(W+Wr)  (5) 
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in  which  u  =  {Ux  uy  uz}T  =  {u  v  w}T  •  Equation  (1)  can  be  solved  for  dis¬ 
placement  or  deflection  of  a  body  subject  to  a  load  described  by  f(x,  y,  z )  and  appropriate 
boundary  conditions.  Displacement  at  specified  points  or  nodes  in  the  domain  of  interest 
are  the  unknown  quantities  or  degree(s)  of  freedom  (dof). 

The  desired  solution,  u(x,y,z),  is  expressed  in  terms  of  the  displacement  vector 
at  each  nodal  point  in  the  domain  and  the  interpolation  functions  within  each  element.  If 
we  define  U  as  the  vector  of  nodal  displacements  with  the  three  components  for  each  node 
grouped  together, 

u(x,y,z)  =  N{x,y,z)U  (6) 

where 

Hi  0  0  H-2  0  0  .  H8  0  0 

N=  0  Hx  0  0  H2  0  .  0  ff8  0  (7) 

0  0  ffj  0  0  if2  .  0  0  h8 

for  a  three-dimensional  hexahedral  element  where  Hi  is  the  three-dimensional  cardinal 
basis  function  corresponding  to  node  i.  Similarly,  the  strain  matrix  [ B ]  is  defined  such  that 


e(x,y,z)  =  B(x,y,z)U  (8) 
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V,  z)  =  [D]B(x,  y,  z)U.  (10) 
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The  dimensions  of  these  elemental  matrices  and  vectors  for  tri-linear  three-dimensional 
hexahedral  elements  are  [U]  =  24  x  1  (8  nodes,  3  degrees  of  freedom  per  node),  [N]  = 
3  x  24,  [B\  =  6  x  24,  and  [D]  =  6  x  6.  U  can  then  be  post-processed  to  calculate  the  stress 
and  strain  fields  as  necessary.  Figure  1  shows  a  notional  hexahedral  finite  element  of  the 
sort  used  here. 


1  -1 


Figure  1 :  Canonical  hexahedral  finite  element 


B.  DISCONTINUOUS  GALERKIN  FORMULATION 

This  section  will  detail  the  formulation  of  nodal  DG  solid  (three-dimensional)  and 
plate  elements  to  be  used  to  solve  for  the  displacement  field  in  a  given  linearly  elastic 
structure.  The  plate  element  to  be  developed  is  based  on  Reissner-Mindlin  theory  but  all 
dof  are  displacements. 

1.  Solid  Element 

Liu,  Wheeler,  and  Dawson  [21]  proposed  and  implemented  a  nodal  DG  formulation 
that  can  be  simply  coupled  with  existing  codes  for  continuous  models  and  can  be  switched 
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between  three  different  specific  methods  via  the  selection  of  a  single  scalar  parameter,  dDG- 
Their  formulation  is  the  starting  point  for  this  work  and  is  summarized  below. 

The  domain,  f)  C  R3,  has  a  boundary,  dQ,  comprised  of  non-intersecting  Dirichlet, 
Tu,  and  Neumann,  T,,  boundaries  upon  which  displacement,  u  e  Hl(Yu),  and  surface 
traction,  t  e  L2(Tt)  are  specified,  x  =  {E\,  E2, . . . ,  EN}  is  a  non-degenerate  discretization 
of  Q;  in  this  work,  E:J  are  hexahedra.  S  =  St  +  Tn  +  T,  is  the  set  of  faces  of  x  where  Si 
are  interior  faces. 

Following  standard  weighted  residual  methodology,  Equation  (1)  is  multiplied  by  a 
test  function  and  integrated  over  the  domain  with  the  intent  of  finding  a  solution  that  leaves 
no  residual  from  that  integral.  Since  this  is  a  Galerkin  formulation,  the  test  function  to  be 
used  is  the  interpolation  function  used  in  the  discretization  of  the  domain;  as  a  Discon¬ 
tinuous  Galerkin  formulation  the  test  function  need  only  be  defined  within  each  element. 
Symbolically, 

f  a(u)  :  S/vdV  —  f  (an)  ■  vdS  —  f  f  ■  vdV.  (11) 

J  E  JdE  J  E 

Taking  advantage  of  the  fact  that 

a(u )  :  Vn  =  a(u)  :  VvT  =  a(u )  :  e(v)  (12) 


the  elemental  equation  in  question  is 


/  a(u)  :  e(v)dV  —  /  (an)  ■  vdS  —If'  vdV. 

IE  JdE  Je 


(13) 


Or,  summed  over  all  elements, 


/  a(y>  ■  £(vW  -  [  (an)  ■ vdS  =  I  f  ■  vdV. 

E&x^E  dE£SddE  EexdE 


(14) 


In  a  continuous  formulation  the  second  term  of  the  last  equation  would  disappear  on  inter- 
elemental  boundaries  and  only  exist  on  the  domain  boundary;  for  DG,  however,  that  is  not 
the  case.  Define  jump  and  average  functions  across  such  an  interior  boundary,  between  two 


10 


elements  arbitrarily  labeled  as  “Left”  and  “Right”  as 


M  =WL~  WR 
1 

W  =  ~{wL  +  wR). 

Combining  those  definitions  with  the  identitiy 

[M  =  +  [<P\W} 


(15) 

(16) 


(17) 


and  the  assumption  that  traction  across  interfaces  is  continuous,  results  in 


E 

dE&{S-ru) ' 


'dE 


{a(w)ns}  •  [v]dS 


Y  [  f-vdV+  Y  [  t-vdS.  (18) 

E&x  E  dE&  Tt  dE 


Liu  et  al.  [21]  add  face  integrals  fgE{a(v)ns}  ■  [u]dS  and  Vp  j.)E [u]  •  [v]dS,  both  of 
which  disappear  for  an  exact  solution,  to  control  symmetry  and  provide  stabilization;  they 
generate  the  following  bilinear  and  linear  forms: 


cr(u)  :  e(v) dV  —  y] 

aEe(Si+r„) 


+  Y,  ®dg 

dE&{Si+Tu) 


{a(u)ns}  ■  [ujdS1 


JdE 

{cr(v)ns}  ■  [rtjdS' 


+  E 

ase(Si+ru) 


SPG 

Isl 


IdE 


u\  •  [u]dS  (19) 
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L(v)  =  £  /  /  •  vdG  +  ^  f  t  ■  vdS  +  ^  f  {a(v)ns}  •  udS 
Eex  E  9E£  rt^dE  dE&Tu  J dE 

+  Y''  f  u  ■  [taels'  (20) 

aEeru  Is!  ^dE 

where  G  is  the  shear  modulus,  <5P  is  a  scalar  penalty  parameter,  \s\  is  the  square  root  of 
the  area  of  the  element’s  face,  and  9Dq  indicates  the  DG  method  in  use.  For  non-zero  5P, 
if  Odg  =  —  1  the  method  is  the  Symmetric  Interior  Penalty  Galerkin  (SIPG),  which  also 
corresponds  to  the  Local  Discontinuous  Galerkin  (LDG)  method  of  Cockburn  and  Shu  with 
f3  =  0  [30];  if  0 do  =  +1  it  is  Non-symmetric  Interior  Penalty  Galerkin  (NIPG);  if  0DG  =  0 
it  is  Incomplete  Interior  Penalty  Galerkin  (IIPG).  If  6dg  =  0  and  dp  =  0  the  method  is  that 
of  Oden,  Babuska,  and  Baumann  (OBB).  The  problem  statement  is  now:  find  u  G  V  such 
that 

a(u,v)  =  L(v)  VveV  (21) 

H\x)  =  {v  e  L'\n)  :  v\Ej  G  H1{Ej)^Ej  Gy};  V  =  {v  G  H\X)}-  (22) 

Equations  (19)  and  (20)  can  be  converted  into  matrix-vector  form.  The  integrand  of 
the  first  term  of  (19)  is: 

cr(u)  :  e(y)  =  DB(x,y,  z)U  :  B(x,y,z )  =  B1  (x,y,  z)DB(x,y,  z)U.  (23) 
Since  U  is  independent  of  (x,  y,  z),  it  can  be  taken  outside  the  integral 

[  a(u)  :  e(v)dV  =  I  BTDBdVU  =  KVU  (24) 

J  E  J  E 

so  that  the  matrix,  Kv  resulting  from  the  volume  integral  is  multiplied  by  the  displacement 
vector.  This  Kv  is  the  identical  to  the  elemental  stiffness  matrix  common  to  continuous 
Galerkin  formulations. 

In  a  similar  fashion,  the  surface  integrals  of  Equation  (19)  can  also  be  expressed 
as  matrix- vector  products.  Once  expanded  the  vector  components  will  be  comprised  of 
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various  products  of  the  unknown  vector,  uL  or  uR  and  the  test  function,  vL  or  vR.  The 
interface  stiffness  matrices  will  be  referred  to  as  KlMN  where  the  superscript,  i  refers  to  the 
integrals  in  the  order  they  appear  in  Equation  (19);  the  subscripts  will  take  the  values  L  and 
R,  referring  to  the  left  and  right  sides  of  the  interface,  respectively.  The  subscript  M  will 
correspond  to  the  v  component  and  the  N  will  correspond  to  the  u  component. 

For  example,  the  expansion  of  the  first  surface  integral  into  the  K\jn  components 
is 


[ {cr(u)ns}  ■  [vJdS'  =  -  f  +  v(ur)  ns)  •  (■ vL  -  vR)dS 

Js  Js  z 

=  f  (cr(uL)ns)  ■  vLdS  +  ^  f  (a (uL)ns)  ■  vRdS 
J  S  J  s 

~  \  [  (a(uR)nS)  •  vLdS  +  ^  f  {a (uR)ns)  ■  vRdS.  (25) 
J  S  s 

Expressing  the  normal  vector  in  terms  of  a  matrix  of  the  direction  cosines, 


As  = 


nx  0  0  riy  0  nz 


0  ny  0  nx  nz  0 

0  0  nz  0  riy  nx 


(26) 


the  terms  of  Equation  (25)  can  be  expressed  as 


2  [WuL)ns)  •  vLdS  =~\f  NTLAsDBLdSUL  =  K\lUl 
J  S  J  s 


(27) 


2  J(v(uL)ns )  •  vRdS  =  NTRAsDBLdSUL  =  KlRLUL  (28) 

-  ^  J  {°{uR)n* )  ■  vLdS  =  ~\Js  NTLAsDBRdSUR  =  K\rUr  (29) 

J  [  (■ <j(uR)ns )  ■  vRdS  =  \  [  NTRAsDBRdSUR  =  KlRRUR.  (30) 
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Conveniently  the  interface  stiffness  matrices  that  result  from  the  third  term  of  (19)  are 
simple  re-arrangements  of  those  calculated  from  the  second  term.  That  is,  K\h  =  Kjj', 
KRR  =  I\RRr ,  K\r  =  K‘rjT ,  and  KR[  =  K\R .  The  final  term  of  (19)  is  referred  to  as 
the  interface  penalty  stiffness,  K\jN  and  its  components  are  calculated  as  follows: 


SpG 

Isl 


u]  ■  [u]dS' 


(uL  -  uR)[yL  -  vR)dS 


(uLvL  -  uLvR  -  uRvL  +  uRvR)dS  (31) 


|s|  Js 

[  uLvLdS  =  ^  [  NlNLdSUL  =  K\lUl  (32) 

M  Js  |s|  Js 

-8jfr  I  uLvRdS  =  [  NTRNLdSUL  =  KrlUl  (33) 

lsl  Js  lsl  Js 

f  u^dS  =  ~T T  f  NTLNRdSUR  =  K3lrUr  (34) 

\S\  Js  M  Js 

^  [  uRvRdS  =Yr  f  NTRNRdSUR  =  KrrUr.  (35) 

lsl  Js  lsl  Js 


The  terms  of  L(v)  are  vectors  resulting  from  integrating  the  body  forces  and  given 
boundary  conditions  (both  displacement  and  traction)  over  the  volume  and  surfaces. 


/  /  •  vdV  —  /  N1  fdV  =  Fb 

'  E  J  E 


where  /  =  (fx,  fy,  fz)T  for  the  element  in  question. 


/  t  ■  vdS  —  /  NrtdS  =  Ft 

'dE  JdE 


/  {o(v)ns}  ■  udS  =  /  B1  D1  (A8)1  udS  =  F'u 

I  dE  JdE 


-y— p  f  U  ■  [V]dS  = 
lsl  JdE  lsl  JdE 

\T  n,  —  „.,\T 


NTudS  =  FP 


where  t  =  (tx,  ty,  tz)T  and  u  =  (u,  v,  w)T  for  the  boundary  face  in  question. 


(36) 


(37) 

(38) 

(39) 
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To  ensure  each  interface  integral  is  calculated  once  and  only  once,  as  the  loop 
through  the  elements  progresses  surface  integrals  are  calculated  only  for  those  faces  corre¬ 
sponding  to  the  positive  canonical  directions.  This  convention  relies  on  element  numbering 
also  proceeding  in  the  positive  directions,  or  element  e’s  +r  neighbor  is  an  element  num¬ 
bered  greater  than  e. 

One  of  the  advantages  of  discontinuous  Galerkin  formulations  is  the  independence 
of  the  elements  with  the  exception  of  the  numerical  fluxes  between  them.  If  this  were 
a  time-dependent  problem  or  if  there  were  some  other  source  of  data  about  neighboring 
elements,  the  various  terms  of  a(u,  v)  could  be  computed  for  each  element  using  the  dis¬ 
placement  data  from  the  previous  time  step  as  the  Ur  terms  for  its  neighbors.  In  a  static 
treatment,  the  various  KlMN  matrices  will  be  assembled  into  a  total  system  stiffness  matrix; 
likewise,  the  various  Fl  vectors  will  be  assembled  into  a  total  system  force  vector  and  a 
global  KU  =  F  and  solved  simultaneously. 

Assembly  of  the  global  [K]  matrix  highlights  the  most  obvious  difference  between 
CG  and  DG  methods:  the  relationships  between  the  elements.  Figures  2  and  3  illustrate 
the  relationships  between  neighboring  elements  in  the  two  methods.  The  interface  between 
the  DG  elements  may  have  the  same  geometric  location  initially,  but  the  faces  are  treated 
as  separate  entities;  each  element’s  dofs  are  wholly  contained  in  that  element.  The  CG 
elements  actually  share  the  face/edge  and  have  common  nodes  and  dofs  there.  Figures  4 
and  5  show  the  connectivity  and  sparsity  patterns  for  the  assembled  global  stiffness  matrices 
of  the  same  four  element  by  two  element  by  one  element  discretization  generated  by  the 
two  methods.  Cursory  examination  of  Figure  5  reveals  the  elemental  connectivity:  the  full 
blocks  along  the  diagonal  are  the  eight  elemental  stiffness  matrices;  the  sparser  off-diagonal 
blocks  represent  the  interfaces.  It  is  apparent,  therefore,  that  element  one’s  neighbors  are 
element  two  and  element  five.  The  connectivity  of  Figure  4  is  related  to  the  individual 
nodes  and  dofs  as  opposed  to  the  elemental  connectivity  of  Figure  5.  These  two  figures 
also  illustrate  one  of  the  significant  disadvantages  of  DG  methods  -  a  tremendous  growth 
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in  the  scope  of  the  numerical  problem;  these  two  matrices  can  be  used  to  solve  the  same 
physical  problem,  yet  the  CG  version  is  only  half  as  large. 
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Figure  2:  Connectivity  between  adjacent  CG  elements.  After  [31] 
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Figure  3:  Connectivity  between  adjacent  DG  elements.  After  [31] 
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Figure  4:  Sparsity  pattern  of  CG  global  stiffness  matrix  for  a  4x2x1  element  structure 
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DG  4x2x1 
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Figure  5:  Sparsity  pattern  of  DG  global  stiffness  matrix  for  a  4x2x1  element  structure 
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a.  Validation 

A  pair  of  beam  problems  were  used  to  validate  code  derived  from  this 
formulation:  a  cantilever  subjected  to  a  concentrated  force  at  the  free  end  and  a  simply- 
supported  beam  subjected  to  a  concentrated  force  at  mid-span.  Both  beams  are  8m  long 
and  have  unit  cross-sectional  area  and  the  same  isotropic  material  properties.  Errors  relative 
to  maximum  deflections  predicted  by  Euler  Beam  theory  were  calculated  for  various  dis¬ 
cretizations  and  plotted  in  Figures  6-7,  both  of  which  demonstrate  quadratic  convergence 
rates  (the  predicted  rate  for  linear  elements)  for  all  three  methods  (NIPG,  SIPG,  IIPG). 


Figure  6:  Convergence  of  three-dimensional  DG  formulation  vs.  discretization  for  can¬ 
tilever  beam  deflection 


2.  Plate  Element 

Adapting  the  above  formulation  to  different  types  of  elements  is  a  matter  of  using 
appropriate  [N],  [B],  and  77  matrices.  Kwon  and  Bang  [32]  developed  a  displacement 
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Figure  7:  Convergence  of  three-dimensional  DG  formulation  vs.  discretization  for  simply- 
supported  beam  deflection 


only  plate  model  following  Mindlin-Reissner  plate  theory  with  the  following  definitions. 

i  I  {ttl  V\  U\-\-n  Vi-^n  W\  U2  V2  U2+n  C 2 .  ;i  W2  ■  ■  •  tin  t)n  ll2n  ^2 n  IV  n\ 

(40) 

where  n  is  the  number  of  nodes  on  the  bottom  of  the  plate,  node  i  +  n  is  taken  to  be  above 
node  i  and  transverse  deflection  of  top  and  bottom  are  taken  to  be  equal,  or  wr  =  wi+n, 
eliminating  the  need  for  wi>n. 


{e&}  -  {ex  ey  q fxy}T  -  [Bb\{U} 

(41) 

Wb}  =  {<JX  CTy  TXy]T  =  [Db]{eb} 

(42) 

{e.}=  {lyz  1XZ}T  =[BS]{U} 

(43) 

{<**}=  {'Tyz  rxz}T  =  [77s]{eJ 

(44) 
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where  the  subscripts  b  and  s  indicate  bending  and  shear,  respectively.  The  volume  integral 
that  is  used  as  the  stiffness  matrix  then  takes  the  form 

[K\=  [  [Bbf[Db][Bb] dft  +  [  [Ba]T[D8][Ba]dn  (45) 


[B 

d=  [[Bn]  [Bb2] 

[Bb3] 

[ Bu ]] 

(46) 

1 

B  | 

0 

H2d~§ 

0 

0 

[Bu]  = 

0 

0 

h2^± 

z  dy 

0 

(47) 

«S  | 

Hi9~§ 

U  dNi 
y 

u  dNi 

0 

[Bs]=  [[Bal\  [Bs2]  [Bs3\  [Bs4\]  (48) 


[Bsi\  = 


0 


0 

OX 


o  Nffi 


n  n  dH2  9Ni 

yj  lint-)  <~v 

L  OZ  OX 


In  the  bending  and  shear  strain  matrices,  [Bb\  and  |  Bs\ ,  II,  (x.  y)  are  the  two  dimensional 
nodal  interpolation  functions  in  the  planar  directions  and  Ni(z)  are  the  one  dimensional 
transverse  nodal  interpolation  functions.  The  shear  term  must  be  under-integrated  numeri¬ 
cally  to  prevent  shear  locking  for  thin  elements.  Further  discussion  of  the  numerical  inte¬ 
gration  schemes  used  can  be  found  in  Appendix  B.  An  approach  used  in  other  DG  formu¬ 
lations  of  the  plate  bending  problem  is  to  use  a  lower  order  interpolant  for  shear  terms  than 
that  used  for  bending  terms  [15]. 


a.  Validation 

A  clamped  square  plate  subject  to  a  concentrated  load  applied  at  its  center 
was  used  to  validate  the  plate  element  model.  The  plate  is  .3048  m  x  .3048  m  x  .00635 
m  (12  in  x  12  in  x  1/4  in)  and  isotropic.  Theoretical  values  were  calculated  according 
to  Timoshenko  [33].  All  three  methods  demonstrate  quadratic  convergence  when  linear 
elements  are  used,  as  shown  in  Figure  8. 
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Figure  8:  Convergence  of  center  deflection  of  a  clamped  plate  comprised  of  plate  elements 

A  second  validation  was  conducted  by  comparing  dynamic  continuous  and 
discontinuous  models  of  the  same  square  plate.  Both  models  consist  of  144  square  elements 
such  that  the  characteristic  length  of  each  element  is  four  times  its  thickness.  Clamped 
boundary  conditions  and  zero  deflection  and  velocity  initial  conditions  were  applied  and 
both  models  have  lumped  (diagonal)  mass  matrices.  Figures  9  and  10  show  excellent  agree¬ 
ment  between  the  two  models  for  the  transverse  displacement  and  velocity  of  the  center  of 
the  plate  subjected  to  a  constant  concentrated  force  applied  at  its  center.  Appendix  A  details 
the  time  integration  algorithms  used. 

One  further  validation  compares  the  present  formulation  with  central  deflec¬ 
tion  of  uniformly  loaded  composite  plates  as  described  by  Lok  and  Cheng  [34].  A  square 
discretization  with  progressively  more  elements  per  side  in  plate  and  a  single  thickness 
element  were  modeled. 
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Figure  9:  Comparison  of  displacement  calculated  for  Continuous  Galerkin  and  Discontin¬ 
uous  Galerkin  models  of  a  center  loaded  clamped  plate 


Figure  10:  Comparison  of  velocity  calculated  for  Continuous  Galerkin  and  Discontinuous 
Galerkin  models  of  a  center  loaded  clamped  plate 
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x  10“3 


Figure  11:  Convergence  of  central  deflection  of  uniformly  loaded  clamped  laminated  plate, 
comparison  with  [34] 
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3.  Stabilization  and  Penalty  Parameters 

One  challenge  in  using  these  discontinuous  Galerkin  elements  is  selection  of  an 
appropriate  penalty  parameter  for  the  last  two  terms  of  Equations  (19)  and  (20).  As  5P 
approaches  infinity  the  methods  return  to  their  continuous  roots,  so  selection  of  too  large 
a  penalty  is  a  waste  of  computing  resources.  The  penalty  must  also  be  large  enough  to 
guarantee  existence  and  uniqueness  of  the  solution  [35].  Additionally,  the  nature  of  the 
penalty  is  often  described  as  a  function  of  the  local  element  order  or  size  with  little  other 
clarification.  The  nature  of  penalizing  the  jump  as  a  stabilization  mechanism  is  discussed 
by  Brezzi  et  al.  [9]  for  elliptic  DG  formulations  in  general.  Others  have  computed  lower 
bounds  on  penalties  for  various  formulations,  [36],  [37],  [38]  among  others. 

In  Liu’s  formulation  [21],  the  penalty  term  is  a  surface  integral  multiplied  by  a 
parameter  to  be  determined  and  the  material  shear  modulus  divided  by  the  square  root  of  the 
area  over  which  the  integral  is  calculated.  Another  intriguing  approach  is  that  of  Ainsworth 
and  Rankin  [39]  in  which  they  compute  a  lower  bound  on  the  penalty  parameter  that  is  a 
function  of  the  method  selector  ( 9dg  above)  and  the  maximum  eigenvalue  of  the  elemental 
stiffness  matrix;  that  value  is  also  divided  by  a  term  analogous  to  the  square  root  of  the  area 
of  the  integral.  In  view  of  this  approach,  the  presence  of  the  shear  modulus  in  Liu’s  penalty 
term  serves  as  a  scaling  factor  to  keep  the  penalty  in  the  same  numerical  neighborhood 
as  that  of  the  volume  integral,  Kv  or  elemental  stiffness  matrix.  Anticipating  using  this 
formulation  for  a  composite  material  with  potentially  vastly  different  shear  moduli  between 
elements,  we  will  combine  these  two  approaches,  replacing  SpG  with  8P-  ma x(p(Kv)).  As 
long  as  5P  >  (1  —  dDG)2,  a  unique  solution  will  exist  [39]. 
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III.  COMPOSITE  STRUCTURES 


Composite  materials  can  provide  designers  with  optimal  combinations  of  strength, 
weight,  flexibility,  and  other  physical  characteristics  for  their  application.  Effective  design, 
however,  requires  a  thorough  understanding  of  material  behavior  in  the  expected  operat¬ 
ing  environment-data  which  can  be  cumbersome  to  obtain,  making  the  development  of 
effective  simulations  key. 

The  strength  of  a  chain  is  determined  by  that  of  its  weakest  link,  but  the  utility  of 
composite  materials  is  the  improvement  each  constituent  element  brings  to  the  desired  ma¬ 
terial  properties  of  the  whole.  Like  alloys,  particulate  composites  generally  demonstrate 
better  properties  than  the  homogeneous  matrix  material  by  virtue  of  the  reinforcement  pro¬ 
vided  by  specifically  chosen  additives. 

A.  MULTI-SCALE  MODEL 

Analysis  of  structures  consisting  of  composite  materials,  and  of  the  materials  them¬ 
selves,  require  multiple  scales  of  analysis  to  be  used  in  conjunction  with  each  other.  That 
is,  the  analysis  can  proceed  from  the  micro-  (atomic  or  molecular)  level  to  the  meso-  (ma¬ 
terial)  level  and  on  to  the  macro-  (structural)  level  and  back  down  again  as  necessary  with 
smearing  or  decomposition  of  properties  and  loads  as  dictated  by  the  analysis  required.  The 
properties  of  a  laminated  fibrous  composite,  for  example,  are  a  function  of  the  properties  of 
both  fiber  and  matrix  as  well  as  the  weave  pattern  of  individual  laminae  and  lay-up  pattern 
of  the  assembled  laminate.  Similarly,  failure  of  the  same  laminate  can  result  from  sepa¬ 
ration  between  laminae,  separation  of  fibers  from  matrix,  or  failure  of  individual  fibers. 
An  illustration  of  this  process  is  shown  in  Figure  12.  Kwon  [40]  presents  an  extensive 
discussion  of  these  cycles. 

This  study’s  focus  on  structural  responses  of  existing  composites  allows  the  ac¬ 
ceptance  of  micro-level  analysis  already  conducted  in  the  design  and  construction  of  the 
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Figure  12:  Multi-scale  analysis  cycle  for  a  fibrous  composite.  From  [40] 


composite  materials  that  will  be  assembled  to  create  the  sandwich  plates  of  interest.  That 
is,  we  will  not  deal  with  the  properties  of  the  fibers  that  are  woven  into  E-glass,  but  will 
instead  treat  the  manufacturer’s  given  orthotropic  properties  of  a  ply  as  known.  Specific 
material  properties  used  are  listed  in  Appendix  B. 

B.  SANDWICH  COMPOSITES 

This  work  will  attempt  to  model  sandwich  composites  for  plate  and  shell  applica¬ 
tions.  These  materials  are  comprised  of  low  density  cores  that  are  relatively  stiff  trans¬ 
versely  and  skin  layers  that  provide  in-plane  strength  to  the  structure.  Because  the  core 
material  is  used  to  provide  transverse  stiffness,  it  is  tempting  to  model  it  using  full  three- 
dimensional  solid  finite  elements;  however,  the  aspect  ratio  of  the  plate  structure  and  its 
included  elements  makes  solutions  generated  using  a  plate  element  for  all  layers  more  ac¬ 
curate.  That  is,  while  the  core  is  generally  the  thickest  component  in  the  sandwich,  it  is 
still  thin  relative  to  its  planar  cross-section,  giving  it  a  sub-optimal  aspect  ratio  for  solution 
with  three-dimensional  elements.  This  model  of  a  composite  material  will  consist  of  an  as- 
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semblage  of  elements  that  are  each  homogeneously  comprised  of  the  constitutent  materials 
and  of  the  type  of  element  described.  That  is,  the  elemental  models  developed  in  Chapter 
II  will  be  used  with  [D]  matrices  describing  the  material  properties  of  constituent  materials 
inserted  appropriately.  Care  must  be  used  to  ensure  the  material  property  matrix  contains 
not  only  the  correct  individual  properties  for  each  element,  but  also  that  it  is  of  the  cor¬ 
rect  form.  Modeling  an  orthotropic  material  with  a  [D]  matrix  calculated  in  isotropic  form 
may  result  in  significant  loss  of  accuracy.  Figure  13  illustrates  the  two  different  transverse 
lay-ups  used  to  model  sandwich  plates  in  this  work;  the  three  layer  model  includes  only 
the  core  and  two  skin  faces;  the  five  layer  or  “with  resin”  model  includes  a  relatively  thin 
layer  with  material  properties  similar  to  common  adhesives  used  in  assembling  sandwich 
composites.  The  two  models  vary  slightly  in  thickness  as  well  as  in  overall  stiffness  due  to 
the  inclusion  of  the  extra  layers. 


Three  layer  plate  model 


skin 


core 


skin 


Five  layer  plate  model 


Figure  13:  Three  and  Five  Layer  Sandwich  Plate  models 

Schmit  and  Monf orton  [41]  developed  a  discrete  element  method  to  predict  the 
static  deflection  of  sandwich  plates  and  shells  with  laminated  faces  under  a  variety  of 
boundary  conditions.  Kanematsu  and  Hirano  [42]  expanded  on  that  work  to  examine  both 
bending  and  vibration  of  sandwich  plates;  their  work  also  included  experimental  validation. 
The  deflection  of  a  50  inch  square  clamped  plate  with  a  one  inch  thick  core  of  aluminum 
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honeycomb  faced  by  two  0.015in  thick  aluminum  skins  under  uniform  pressure  was  calcu¬ 
lated  by  both  papers  and  is  used  here  to  validate  the  DG  structural  model  developed  in  the 
previous  chapter.  Using  one  thickness  element  for  each  face  and  one  for  the  core,  maxi¬ 
mum  static  deflection  was  calculated  for  a  range  of  square  planar  discretizations;  they  are 
plotted  in  Figure  14.  Good  agreement  was  reached  with  as  few  as  four  elements  in  each 
direction,  and  refinement  further  than  twelve  elements  in  each  direction  was  shown  to  be 
unnecessary.  Both  sets  of  authors  neglected  in-plane  bending  of  the  core  material,  which 
the  present  formulation  does  not,  therefore  the  current  model  is  sitffer  and  returns  a  slightly 
lesser  static  deflection. 

Relative  convergence  of  both  maximum  deflection  and  maximum  bending  stress 
in  a  clamped,  square,  five-layer  sandwich  plate  subject  to  a  concentrated  center  force  was 
examined  by  modeling  a  quarter  of  the  plate,  taking  advantage  of  symmetry  to  achieve 
finer  discretizations  without  incurring  excessive  computational  cost.  For  this  study,  a  range 
of  twelve  to  thirty  elements  per  side  of  the  quarter-plate  was  modeled.  Deflection  is  the 
primary  variable  and  bending  stresses  are  post-processed  quantities.  Both  skin  and  core 
stresses  were  calculated,  but  as  there  was  minimal  difference  between  the  two,  core  stresses 
are  omitted  from  the  plot  for  clarity.  Convergence  was  calculated  relative  to  the  finest 
discretization  calculated,  thirty  elements  per  quarter-plate  side  (dx=0.0075  m)  and  is  shown 
in  Figure  15.  Both  quantities  converge  at  better  than  quadratic  rates,  the  predicted  rate  for 
linear  elements. 


30 


0.1 


0.08 


—  0.06 
L 

T3 


0.02 


calculated  max  deflection 
predicted  max  deflection 


5  10  15  20 

number  of  elements  per  side 


25 


Figure  14:  Convergence  of  DG  Sandwich  Plate  max  deflection  (due  to  uniform  load)  to 
predicted,  comparison  with  [41,  42] 


Figure  15:  Relative  convergence  of  DG  Sandwich  Plate  max  deflection  and  skin  stress  (due 
to  concentrated  load) 
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c. 


FAILURE  MODES  AND  IDENTIFICATION 


Failure  of  sandwich  structures  is  a  function  of  the  constitutent  materials,  the  ge¬ 
ometry  of  the  structure,  and  the  nature  of  the  loading.  Common  failure  modes  of  such 
structures  include  debonding,  delamination,  core  crushing,  skin  wrinkling,  and  general 
buckling.  Debonding  is  the  separation  of  the  skin  material  from  the  core  and  delamination 
usually  refers  to  the  separation  of  layers  within  the  skin  material.  In  this  work,  debonding 
will  be  the  primary  failure  mode  examined.  No  attempt  is  made  here  to  develop  failure 
criteria;  we  will  instead  attempt  to  develop  an  appropriate  method  to  reflect  failure  within 
a  structural  model. 

Previous  work  [43]  concluded  that  modeling  an  independent  layer  representing  the 
adhesive  between  laminae  of  composites  was  necessary  to  observe  the  delamination  failure 
mode.  In  this  work  we  will  examine  whether  such  a  layer  can  be  omitted  when  Discontin¬ 
uous  Galerkin  (DG)  techniques  are  applied  in  assembling  the  structure. 

One-quarter  of  a  twenty-four  by  twenty-four  element  square  clamped  plate  with  the 
same  aluminum  skins  and  honeycomb  core  as  was  used  in  the  previous  section  was  the 
basis  for  this  examination.  This  plate  model  is  450  mm  x  450  mm,  the  core  is  10  mm 
thick,  each  skin  is  0.375  mm  thick,  and  the  load  is  a  concentrated  force  of  1000  N  applied 
to  the  center  of  the  plate.  The  model  was  assembled  as  described  in  Chapter  II  and  the 
global  displacement  vector  was  calculated.  The  degrees  of  freedom  of  interest,  those  on 
the  interface  between  the  bottom  of  the  core  and  the  top  of  the  lower  skin,  were  extracted 
and  post-processed  to  calculate  the  bending  stress  vector  at  each  point  on  that  interface. 
In  order  to  display  the  calculated  data,  an  interpolation  function  describing  resulting  stress 
values  of  interest,  crx,  was  generated  using  MATLAB’s  TriScatteredlnterp  function 
and  then  applied  over  a  grid  of  the  same  dimension  as  the  original  discretization  and  plotted 
using  contourf. 

Figures  16  and  17  show  the  resulting  normal  stress  in  the  x  direction  on  the  top 
of  the  lower  skin  (the  side  facing  the  core)  and  on  the  bottom  of  the  core  (the  side  facing 
the  lower  skin)  for  a  three  layer  model.  The  upper  right-hand  comer  of  these  quarter-plate 
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plots  corresponds  to  the  center  of  the  whole  plate  and  is  the  point  of  application  of  the 
concentrated  load.  As  expected,  the  peak  stress  values  for  both  components  are  located  at 
the  center  of  the  plate,  and  the  stiffer  skin  is  taking  a  significantly  larger  portion  of  this  in¬ 
plane  load.  The  effects  of  the  clamped  boundary  conditions  can  also  be  seen  at  the  edges  of 
both  graphs.  This  process  was  repeated  for  a  five-layer  model  with  a  lay-up  of:  skin-resin- 
core-resin-skin.  In  this  case  the  resultant  stresses  were  calculated  for  the  top  of  the  lower 
skin  and  bottom  of  the  core  as  before  as  well  as  both  the  top  and  bottom  of  the  intervening 
resin  layer.  The  results  are  displayed  in  Figures  18  -  21.  The  stresses  on  the  two  faces 
of  the  resin  layer  are  close  enough  in  magnitude  to  treat  them  as  equal.  Additionally,  the 
stresses  on  the  skin  and  core  are  negligibly  affected  by  the  insertion  of  the  resin  layer  into 
the  model.  Unfortunately,  the  relative  magnitudes  of  the  stresses  do  not  suggest  an  intuitive 
or  convenient  stress-based  failure  criteria  that  could  be  applied  to  the  resin  in  absentia 
based  on  the  calculated  values  on  the  faces  of  the  skin  and  core.  Therefore,  the  inclusion 
of  an  interface  layer  is  needed  to  accurately  model  debonding  of  skin  faces  from  cores  of 
sandwich  composites. 


Figure  16:  Planar  stress  on  the  skin  of  a  clamped  three-layer  sandwich  plate  subject  to 
concentrated  center  force 
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Figure  17:  Planar  stress  on  the  core  of  a  clamped  three-layer  sandwich  plate  subject  to 
concentrated  center  force 


Figure  18:  Planar  stress  on  the  skin  of  a  clamped  five-layer  sandwich  plate  subject  to 
concentrated  center  force 
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Figure  19:  Planar  stress  on  the  core  of  a  clamped  five-layer  sandwich  plate  subject  to 
concentrated  center  force 


Figure  20:  Planar  stress  on  the  bottom  of  a  resin  layer  of  a  clamped  five-layer  sandwich 
plate  subject  to  concentrated  center  force 
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Figure  21:  Planar  stress  on  the  top  of  a  resin  layer  of  a  clamped  five-layer  sandwich  plate 
subject  to  concentrated  center  force 
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D.  FAILURE  MODELING 


This  section  will  compare  and  contrast  three  proposed  methods  of  modeling  debond¬ 
ing  failure  between  the  core  and  the  resin  layer  opposite  an  imposed  concentrated  force  at 
the  center  of  a  sandwich  plate.  Each  of  the  methods  to  follow  will  be  demonstrated  using  a 
five  layer  plate  discretized  into  a  twenty-four  by  twenty-four  element  mesh.  The  calculated 
planar  stress  values  on  the  bottom  of  the  core,  the  top  and  bottom  of  the  lower  resin  layer 
and  the  top  of  the  skin  layers  will  be  displayed  and  discussed.  The  undamaged  model  as 
displayed  in  Figures  18-21  will  be  used  as  a  baseline  reference  case. 

1.  Damage  via  Complete  Disconnection 

Mergheim  et  al.  [22]  introduce  a  scheme  that  combines  DG  methods  with  existing 
interface  methods  for  modeling  failure.  Specifically,  they  propose  re-defining  Equation 
(19)  as  follows: 

a(u,v)  =  '<^2  f  cr(u)  :  e(v) 6V  —  (1  —  a)  f  { a (u) ns }•  [njdS1 

egx  ^ e  dE&(Si+ r„)  ^ dE 

+  ^  #dg(1  -  «)  f  {a(v)ns}  •  [t^dS 

dE&(Si+ r„)  JdE 

+  f  [n]  •  ^(1  —  a)^yy[u\  +  at[u]^  dS.  (50) 

dEe(Si+ru)  J dE  ^  Is!  ' 

where  a  is  a  switching  factor  and  t  is  a  traction  vector  governed  by  a  traction-separation 
law.  In  the  pre-critical  or  undamaged  regime,  a  =  0  and  Equation  (50)  is  identical  to 
Equation  (19);  in  the  post-critical  or  damaged  regime,  a  =  1  and  the  surface  integrals 
representing  connected  interfaces  are  replaced  by  a  traction-separation  law  that  models 
progressive  failure. 

Adapting  this  concept  and  assuming  complete  failure  of  the  interface  between  the 
core  and  resin,  debonding  damage  was  simulated  by  separating  a  four  by  four  array  of 
resin  elements  surrounding  the  center  of  the  plate  from  their  core  element  neighbors.  In 
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the  present  formulation,  all  of  the  area  integrals  of  Equation  (19)  are  calculated  for  both 
interior  and  Dirichlet  exterior  boundaries.  If  the  interface  between  two  elements  is  deemed 
to  have  failed,  those  faces  can  then  be  considered  members  of  Tu  rather  than  Si,  so  the 
K)j  and  KRR  terms  remain  and  the  various  left  and  right  components  {K\R  and  KRL)  are 
simply  deleted  from  the  global  K  matrix. 

All  other  parameters  were  unchanged  from  the  undamaged  case.  The  resulting 
stress  fields  are  displayed  in  Figures  22-25.  The  skin  stresses  are  largely  unchanged  from 
the  undamaged  state.  The  core,  however,  bears  the  brunt  of  this  simulated  debonding;  its 
maximum  stress  value  is  twenty  times  that  of  its  undamaged  version.  This  is  because  this 
simulation  has  effectively  removed  all  constraints  on/supports  to  the  core  in  the  area  of 
greatest  load.  In  the  undamaged  model,  the  static  deflection  at  the  center  of  the  plate  is 
consistent  through  the  thickness;  that  is,  the  transverse  displacement  at  all  center  nodes, 
in  both  skin  layers,  both  resin  layers,  and  the  core,  has  been  the  same.  In  this  example, 
the  static  deflection  is  consistent  from  the  top  of  the  structure  (point  of  load  application) 
down  to  the  bottom  of  the  core;  the  deflection  of  the  lower  resin  and  lower  skin  layers  was 
consistent  within  those  two  layers,  but  markedly  less  than  that  above.  Deflection  curves 
along  the  centerline  of  the  plate  for  the  bottoms  of  the  core,  resin,  and  lower  skin  layers  are 
shown  in  Figure  26.  This  is  because  this  method  of  debonding  the  plate  has  also  eliminated 
any  means  of  transferring  the  load  between  those  layers  within  the  damage  zone.  Complete 
disconnection  of  inter-elemental  interfaces  is  not  the  proper  way  to  model  this  sort  of  dam¬ 
age.  A  more  consistent  approach  may  be  to  disconnect  the  planar  dofs  between  elements, 
but  leave  the  transverse  dofs  connected. 
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Figure  22:  Planar  stress  on  the  skin  of  a  damaged  via  complete  disconnection  clamped 
five-layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  23:  Planar  stress  on  the  core  of  a  damaged  (via  complete  disconnection)  clamped 
five-layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  24:  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (via  complete  dis 
connection)  clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 


Figure  25:  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (via  complete  discon 
nection)  clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  26:  Static  deflection  along  centerline  of  damaged  (via  complete  disconnection) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 
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2.  Damage  via  Partial  Disconnection 

To  implement  a  more  physically  consistent  disconnection  between  core  and  resin 
elements  in  the  debonding  zone,  only  those  rows  and  columns  of  the  interface  sub-matrices 
that  correspond  to  the  planar  (ut  and  vt)  dofs  will  be  removed.  That  is,  the  entries  in  those 
sub-matrices  that  correspond  to  the  transverse  (wy)  dofs  will  be  left  in  place.  Once  again, 
this  removal  of  entries  from  the  global  stiffness  matrix  is  executed  after  its  assembly  -  a 
step  that  can  be  repeated  as  necessary  for  progressive  failure  with  relative  simplicity.  The 
resultant  stresses  can  be  seen  in  Figures  27  -  30.  The  maximum  stress  values  are  consis¬ 
tently  lower  than  the  undamaged  case  in  all  three  materials,  but  more  noteworthy  is  the 
movement  of  the  location  of  maximum  stress  from  the  center  of  the  domain,  directly  below 
the  imposed  load,  to  the  edge  of  the  debonding  zone.  This  behavior  appears  to  be  consistent 
with  a  physical  stress  concentration  on  the  edge  of  a  discontinuity  in  a  structure  and  may 
be  useful  in  modeling  damage  propagation.  Figure  31  shows  that  the  static  deflection  cor¬ 
responding  to  this  failure  model  is  physically  consistent:  the  core  does  not  deflect  beyond 
or  through  the  resin  and  skin  layers  below  it,  but  all  three  layers  show  sharper  deflection 
within  the  debonding  zone  than  in  the  rest  of  the  domain. 

Figures  32  -  34  show  the  stress  profile  generated  using  this  partial  disconnection 
method  with  the  quarter-plate  model  and  a  discretization  of  thirty  elements  per  side.  The 
stress  concentration  effect  observed  in  the  relatively  coarse  meshes  of  Figures  27  -  30  are 
clearly  present  in  the  finer  model  as  well. 
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Figure  27:  Planar  stress  on  the  skin  of  a  damaged  (by  partial  disconnection  clamped)  five 
layer  sandwich  plate  subject  to  concentrated  center  force 


Figure  28:  Planar  stress  on  the  core  of  a  damaged  (by  partial  disconnection  clamped)  five 
layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  29:  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (by  partial  discon¬ 
nection)  clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 


Figure  30:  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (by  partial  disconnection) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  31:  Static  deflection  along  centerline  of  damaged  (via  partial  disconnection)  five 
layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  32:  Planar  stress  on  the  skin  of  a  damaged  (by  partial  disconnection)  clamped  five 
layer  sandwich  plate  subject  to  concentrated  center  force,  fine  view 
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Figure  33:  Planar  stress  on  the  core  of  a  damaged  (by  partial  disconnection)  clamped  five- 
layer  sandwich  plate  subject  to  concentrated  center  force,  fine  view 


Figure  34:  Planar  stress  in  the  resin  layer  of  a  damaged  (by  partial  disconnection)  clamped 
five-layer  sandwich  plate  subject  to  concentrated  center  force,  fine  view 
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3.  Damage  via  Reduced  Moduli 

The  next  method  of  imposing  damage  in  this  model  was  to  leave  all  elements  con¬ 
nected  as  in  the  undamaged  state,  but  to  reduce  the  effectiveness  of  the  resin  elements  in 
the  damage  zone.  The  same  resin  elements  identified  in  the  previous  attempt  remained 
connected  to  their  core  counterparts,  but  their  shear  and  elastic  moduli  were  reduced  to  1% 
of  the  values  used  for  the  rest  of  the  layer-an  arbitrarily  chosen  reduction.  The  resultant 
stresses  are  plotted  in  Figures  35-38.  In  this  model,  the  stress  in  the  core  elements  is  com¬ 
parable  to  that  of  the  undamaged  plate  and  the  dramatically  lower  stress  values  in  the  center 
of  the  resin  layer  is  entirely  attributable  to  the  lower  modulus.  This  method  would  seem  to 
be  similar  to  Mergheim’s  insertion  of  a  traction-separation  law  to  describe  the  progression 
from  damage  initiation  to  complete  separation  [22] . 


Figure  35:  Planar  stress  on  the  skin  of  a  damaged  (by  reduced  modulus)  clamped  five-layer 
sandwich  plate  subject  to  concentrated  center  force 
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Figure  36:  Planar  stress  on  the  core  of  a  damaged  (by  reduced  modulus)  clamped  five-layer 
sandwich  plate  subject  to  concentrated  center  force 


Figure  37 :  Planar  stress  on  the  bottom  of  the  resin  layer  of  a  damaged  (by  reduced  modulus) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 
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Figure  38:  Planar  stress  on  the  top  of  the  resin  layer  of  a  damaged  (by  reduced  modulus) 
clamped  five-layer  sandwich  plate  subject  to  concentrated  center  force 
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E.  SYNTHESIS 


Table  1  contains  a  summary  of  the  above  results;  the  numerical  stress  values  dis¬ 
played  are  of  use  only  in  their  relationship  to  each  other.  These  were  all  calculated  from 
the  twelve  by  twelve  element  quarter-plate,  which  corresponds  to  the  coarsest  discretization 
used  in  the  convergence  calculations  displayed  in  Figure  15. 

Complete  elimination  of  interface  terms  from  the  global  stiffness  matrix  does  not 
correctly  model  physical  constraints  on  the  core  from  remaining  layers  in  the  structure  be¬ 
low  the  section  deemed  to  have  been  delaminated.  Disconnection  of  the  planar  terms  whilst 
retaining  transverse  interface  terms  does  appear  to  correctly  model  expected  physical  be¬ 
havior.  Adjusting  the  physical  properties  of  the  resin  or  interface  layer  remains  a  potential 
tool  for  faithful  modeling  of  traction-separation  laws. 


Undamaged 

Complete 

Disconnect 

Semi- 

Disconnect 

Reduced 

Modulus 

value 

location 

value 

location 

value 

location 

value 

location 

skin 

96.88 

center 

77.89 

center 

41.02 

zone  edge 

97.61 

center 

core 

0.280 

center 

4.389 

center 

0.111 

zone  edge 

0.284 

center 

resin  top 

10.82 

center 

8.639 

center 

4.413 

zone  edge 

4.331 

zone  edge 

resin  bottom 

10.85 

center 

8.724 

center 

4.641 

zone  edge 

4.303 

zone  edge 

Table  1 :  Comparison  of  maximum  planar  stress  values  (in  MPa)  and  locations  for  damaged 
and  undamaged  clamped  sandwich  plates  subject  to  a  concentrated  center  force 
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IV.  FLUID  MODEL 


Following  the  work  of  Olson  and  Bathe  [23],  Fluid- Structure  Interaction  (FSI)  will 
be  analyzed  via  a  velocity  potential  in  the  fluid.  The  transverse  plate  velocity  will  be 
matched  to  the  z  component  of  the  fluid  velocity  as  a  compatibility  condition  between  the 
two  domains.  The  scalar  velocity  potential  follows  the  wave  equation,  so  a  sufficiently  ac¬ 
curate  and  efficient  model  of  that  equation  with  appropriate  boundary  conditions  is  required 
for  the  fluid  portion  of  this  work. 

A  great  many  methods  are  available  to  model  the  fluid  mechanics  necessary  for  this 
work,  but  as  our  primary  interest  is  in  the  structural  side  of  the  fluid- structure  interaction, 
achieving  adequate  accuracy  without  incurring  significant  computational  cost  lead  to  the 
exploration  of  CA  methods. 

A.  VELOCITY  POTENTIAL  AND  WAVE  EQUATION  THEORY 

The  velocity  potential,  q i>,  in  the  fluid  domain  is  defined  as 

v  =  V0  (51) 

where  v  is  the  velocity  of  the  fluid.  The  wave  equation  is 

u  =  c2V2u  (52) 

where  c  is  the  acoustic  speed  of  the  (fluid)  medium.  Coupled  with  appropriate  initial  con¬ 
ditions,  this  well-posed  initial  value  problem  has  been  much  studied  and  discussed  [44].  In 
one-dimension,  the  problem 


Utt  =  C2UXX  —  OO  <  X  <  oo  0  <  t  <  oo 


«(z,0)  =  f(x) 


(53) 
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Ut(x,  0)  =  g(x) 


is  satisfied  by  the  D’Alembert  solution 

1  1  rx+ct 

t)  =  - [f(x  -  ct )  +  f(x  +  ct )]  +  —  /  g{£)d£.  (54) 

"  J x—ct 

The  three-dimensional  extensions  are 


utt  =  c2V2u  (x,  y,  z )  e  R3 


(55) 


u(x,y,z,  0)  =  <p(x,y,z) 
ut(x,y,z,  0)  =  il>(x,  y,  z) 

which  is  satisfied  by 

d 

u(x,y,z,t)  =  tip  +  —  [t<p\  (56) 

where  ip  and  (p  are  the  averages  of  their  respective  initial  conditions  over  the  sphere  of 
radius  ct  centered  at  (x,  y,  z);  specifically, 


f*  7T  p2n 


x,y,z )  = 


4nc2t2 


'o  Jo 


ip(x  +  ct  sin  4 b  cos  9,y  +  ct  sin  4>  sin  9 , 

z  +  ct  cos  4>)  (ct)2  sin  (pd9d(p  (57) 


and 


7T 


<p(x,  y,  z)  =  ,o2  /  /  <p(x  +  ct  sin  4>  cos  9,  y  +  ct  sin  4>  sin  9, 

Anc~t  J 0  J 0 


z  +  ct  cos  (p)(ct)2  sin  <pd9d(p.  (58) 


The  integrals  in  (57)  and  (58)  are  rarely  simple  to  evaluate  analytically. 
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B.  FINITE  ELEMENT  MODEL  OF  THE  WAVE  EQUATION 


Beginning  with  the  wave  equation  as  applied  to  the  velocity  potential, 

2_9  ,  d2(j)  •• 

c  v  * = w  =  * 

multiply  all  terms  by  a  test  function,  w,  and  integrate  over  the  domain  to  get 

f  wcp dO  —  c2  f  wV20 dO  =  0 
Jn  Jn 

integrate  the  second  term  by  parts  and  apply  Green’s  Identity  to  get 


/  w(J) dfl  +  c2  /  Vw  ■  VcidO  =  c2  /  tcV0  •  hdr  (61) 

Jn  Jn  Jr 

recall  the  definition  of  the  velocity  potential  to  transform  the  right-hand  side, 


/  wcpdtt  +  (?  /  VwVcpdQ  =  c2  /  wn-ndT.  (62) 

Jn  Jn  Jr 

Choosing  Galerkin  test  functions,  the  derivatives  of  the  trial  functions  with  respect  to  the 
unknowns,  the  two  volume  integrals  become 


[M,\  =  /  {i/}T{i/}dfi 
Jn 


[K,\  =  /  {SJH}t{VH} <K2 
Jn 


where  { H }  is  the  vector  of  nodal  interpolation  functions.  In  the  usual  Finite  Element  (FE) 
matrix-vector  form  we  get 


[Mf]{(j)}  +  c2[Kf]{(j)}  =  c2  J  wv  ■  ndr. 
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Equations  (62)  and  (65)  hold  for  each  element  of  a  finite  element  domain,  and,  if 
continuity  between  elements  is  enforced,  also  hold  globally.  Thus  the  right-hand  side  is 
only  defined  on  the  boundary  of  the  domain.  At  the  fluid- structure  interface  the  velocity 
compatibility  provides  a  convenient  input  to  this  finite  element  problem.  Specified  Dirichlet 
or  Neumann  boundary  conditions  can  also  be  applied  with  relative  simplicity,  but  for  this 
work  non-reflecting  boundary  conditions  are  most  appropriate  yet  are  not  easily  applied. 

C.  CELLULAR  AUTOMATA  MODEL  OL  THE  WAVE  EQUATION 

Cellular  Automata  (CA)  are  discrete,  rule-based  numerical  methods  that  can  model 
complex  physical  phenomena  with  relative  simplicity.  Generally,  both  space  and  time  are 
treated  discretely  and  the  value  of  the  quantity  in  question  is  limited  to  a  finite  set  of  val¬ 
ues.  As  the  space-time  domain  proceeds  or  grows  the  seemingly  simple  model  converges 
to  the  complex  real-world  behavior.  The  simplicity  of  the  chosen  rules  and  their  imple¬ 
mentation  lowers  the  computational  cost  while  still  achieving  required  accuracy.  CA  rules 
developed  for  modeling  wave  propagation  are  pre-cursors  to  the  lattice  Boltzmann  method 
of  modeling  fluid  flow. 

Following  the  work  of  Chopard  [25,  26],  Kwon  and  Hosoglu  [28]  modeled  the  wave 
equation  in  one-  and  two-dimensions  with  fixed  boundaries  using  the  following  rules: 

+  At)  =  (f>w(t)  ~  -  At)  +  (f>E(t)  (66) 

4>c{t  +  At)  =  ( )  +  <f>E(t)  +  <f>s(t)  +  <pN(t)  -  2 <f>c(t  -  At))/2.  (67) 

The  value  of  0  at  each  interior  grid  point  in  the  domain  of  interest  (0c)  is  updated  according 
to  the  values  at  its  nearest  Von  Neumann  neighbors  (0at,  0s,  4>e,  4>w)  as  shown  in  Figure 
39. 

For  convenience  the  set  of  points  is  divided  into  two  sets,  “black”  and  “white”  (or 
“odd”  and  “even”),  such  that  the  neighbors  of  each  white  point  are  all  black,  and  only  one 
“color”  is  updated  during  each  iteration.  This  model  includes  an  evolution  of  the  variable  0 
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Figure  39:  Center  node  (black)  and  its  Von  Neumann  neighbors  (white) 

from  being  restricted  to  integer  values  (as  in  a  traditional  CA  model)  to  being  real  valued. 
In  its  fully  discrete  form  this  CA  rule  was  developed  to  model  particle  motion;  with  the 
relaxation  that  allows  real-valued  states,  it  also  corresponds  to  the  finite  difference  model 
of  the  wave  equation  on  a  uniform  grid. 

1.  One  Dimension 

The  classic  illustration  of  D’Alembert’s  solution  to  the  one-dimensional  wave  equa¬ 
tion  is  the  perturbation  of  a  string  subject  to  tension.  For  the  moment,  we  shall  apply  fixed 
boundary  conditions  to  the  ends  of  the  “string”  and  focus  our  attention  to  interior  points 
well  away  from  those  ends.  Further  discussion  of  appropriate  boundary  conditions  will 
follow.  Consider  a  string  subject  to  a  Gaussian  perturbation  at  its  midpoint  as  shown  in 
Figure  40, 


Utt  =  C2UXX  —  OO  <  X  <  oo  0  <  t  <  oo 

_  A 

u(x,  0)  =  f(x)  =  e  2 


(68) 
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ut(x,  0)  =  g(x)  =  0. 


The  D’Alembert  solution  to  this  problem  is 


1.  (x  —  ct)2  ( x+ct )2 

u[x,t)  —  ~[e~  9  +e  ® 


(69) 


which  agrees  nicely  with  the  CA  solution  computed  using  the  rule  found  in  Equation  (66) 
as  shown  in  Figure  41. 


Figure  40:  Initial  perturbation  of  an  infinite  string 
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Figure  41:  CA  solution  vs.  D’Alembert’s  solution  to  the  one-dimensional  wave  equation 
in  an  infinite  string 

2.  Boundary  Conditions 

Fixed  boundary  conditions  are  easy  to  implement,  but  are  of  limited  utility  in  mod¬ 
eling  a  potentially  infinite  fluid  domain.  One  possible  approach  is  to  model  a  much  larger 
domain  than  that  of  interest  so  that  the  area  of  interest  is  well  away  from  the  boundary  and, 
as  such,  solutions  within  it  are  unpolluted  by  whatever  boundary  condition  is  imposed.  This 
requires  many  computations  that  will  be  ignored-a  seeming  waste.  Another  approach  is  to 
generate  a  model  of  a  non-reflecting  boundary  such  that  the  wave  in  question  is  unaffected 
by  its  proximity.  This  is  an  active  area  of  research  in  the  finite  element  arena.  In  the  CA 
arena,  Chopard  and  Droz  [25]  suggested  implementing  boundary  conditions  by  generat¬ 
ing  virtual  cells  adjacent  to  the  boundary  cells  as  shown  in  Figure  42.  Fixed  or  specified 
boundary  conditions  do  not  require  a  virtual  neighbor,  but  are  shown  for  completeness. 
Non-reflecting  (or  zero-gradient  or  adiabatic)  boundaries  model  the  behavior  in  the  heart 
of  the  domain  of  interest,  well  away  from  the  influence  of  any  boundary.  The  free  boundary 
condition  corresponds  to  that  of  a  constant  gradient.  In  practice,  this  can  be  implemented 
either  by  generating  the  virtual  neighbor  cells  or  by  applying  the  resulting  rule  to  the  actual 
cells  located  on  the  boundary  in  question. 
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fixed 


Figure  42:  Virtual  cell  values  for  various  boundary  conditions  in  one  dimension.  After  [25] 

Figure  43  shows  the  application  of  several  different  boundary  conditions  to  the 
positive  side  of  the  spatial  domain  of  our  one-dimensional  wave  equation.  For  reference 
the  f(x  —  ct)  portion  of  the  D’Alembert  solution  is  shown  as  an  exact  solution.  The  non¬ 
reflecting  boundary  condition  corresponds  quite  well  with  the  analytic.  All  four  waves  peak 
initially  in  unison;  the  reflecting  wave  (black)  returns  with  equal  amplitude  while  the  free 
wave  (red)  returns  with  an  inverse  amplitude.  The  fluid  domain  to  be  modeled  will  have 
non-reflecting  boundary  conditions  imposed  at  the  arbitrary  edge  of  the  domain  and  a  free 
boundary  used  to  represent  the  air- water  interface. 

3.  Discretization  and  Model  Fidelity 

The  various  waves  modeled  above  by  CA  rules  display  a  coarseness  that  results 
from  updating  the  value  of  each  point  in  space  at  alternating  time  steps.  The  initial  per¬ 
turbation  displayed  thus  far  has  been  a  Gaussian  wave  of  medium  width  that  proved  to 
be  smooth  enough  to  demonstrate  the  desired  characteristics.  Attempts  to  model  a  point 
source  along  the  lines  of  f(x)  =  0,  x  ^  0,  f(x)  —  1,  x  —  0  were  unsuccessful,  begging  the 
question  of  how  smooth  a  function  or  discretization  are  necessary  to  use  CA  to  model  the 
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traveling  wave  vs  time  at  a  common  spatial  point 


Figure  43:  Application  of  various  boundary  conditions  to  CA  calculation  of  one¬ 
dimensional  wave  propagation 


wave  equation.  In  order  to  determine  a  sufficiently  fine  discretization  relative  to  the  sharp- 

/  \  -A 

ness  of  the  initial  perturbation,  f(x)  =  e  ^ ,  a  series  of  progressively  narrower  Gaussians 
as  shown  in  Figure  44.  For  each  initial  condition,  the  CA  solution  to  the  wave  equation  was 
calculated  for  the  same  set  of  discretizations  (progressively  smaller  dx)  and  error  norms  rel¬ 
ative  to  the  D’Alembert  solution  were  calculated  at  the  same  arbitrary  time.  Non-reflecting 
boundary  conditions  were  applied  in  all  cases.  The  resultant  convergence  as  a  function  of 
dx  is  shown  in  Figure  45.  An  error  norm  of  1%  was  chosen  as  the  comparison  point.  The 
largest  dx  value  required  to  achieve  that  level  of  convergence  was  then  plotted  against  the 
width  at  half  maximum  for  each  initial  condition  as  shown  in  Figure  46.  A  linear  estimate 
of  that  data  is  that  dxc  —  |cr.  Applying  this  estimate  to  the  specific  f(x)  used,  4.5  elements 
or  nodes  are  required  to  represent  the  descent  from  peak  value  of  /  to  1%  of  that  peak. 
Therefore,  a  discretization  that  uses  eleven  or  more  nodes  to  represent  both  sides  of  a  peak 
should  be  sufficient. 


59 


1 


0.8 


0.6 


o 


0.4 


0.2 


-10 


-5 


0 

x 


exp(-x2/2) 

exp(-x2) 

exp(-2x2) 

exp(-4x2) 

exp(-8x2) 


10 


Figure  44:  Initial  perturbations  of  varying  widths 
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dx 


Figure  45:  CA  solution  to  Id  wave  equation  convergence  as  a  function  of  dx  for  various 
initial  conditions 
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dx  ...  ,  vs  full-width  at  half-maximum 
critical 


Figure  46:  Critical  discretization  vs.  initial  perturbation  width 
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4.  Convergence 

To  examine  convergence  of  the  present  CA  rule  for  wave  propagation  with  respect 
to  mesh  or  lattice  spacing,  comparison  to  the  steady-state  plane  wave  as  presented  by  Junger 
and  Feit  [45]  was  used.  Consider  a  semi-infinite  fluid-filled  space  with  a  given  uniform 
vibration,  w{t)  =  W e~lult  at  the  z  =  0  boundary  and  a  non-reflecting  boundary  as  z 
oo.  The  steady  state  pressure  in  the  wave  guide  is  p(z,t)  =  pcW e^kz~lU}t^  where  k  = 
The  given  function  was  applied  to  the  z  —  0  nodes  in  a  CA  domain  with  the  initial 
values  everywhere  else  uniformly  zero,  the  CA  rule  was  applied  for  a  number  of  iterations 
corresponding  to  over  three  periods  of  the  steady-state  solution,  and  point-by-point  error 
was  calculated  relative  to  the  analytic  solution  over  the  range  z  G  [0,  \  and  plotted  in 

Figure  47.  The  rate  of  convergence  is  linear,  as  expected  for  a  first-order  method. 
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Figure  47 :  Convergence  of  C A  wave  equation  rule  to  analytic  solution  as  a  function  of  dx 
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5. 


Three  Dimensions 


It  follows  from  Equations  (66)  and  (67)  that  the  three-dimensional  wave  equation 
can  be  modeled  as 


4>c(t  +  At)  =  ~[<t>w{t)  +  +  <t>s(t)  +  <i>N(t)  +  <M*)  +  <M*)  -  3 <j>c{t  -  At)]  (70) 


with  the  subscripts  on  the  new  terms  standing  for  “front”  and  “back,”  respectively.  To  test 
this  supposition,  consider  a  “point”  source  located  in  the  center  of  a  domain  of  interest. 
As  was  discussed  above,  CA  is  not  expected  to  faithfully  model  a  true  point  source,  so 
the  source  under  consideration  is  a  smooth  radial  function  that  has  a  maximum  value  of 
1  at  the  origin  of  the  domain  and  is  zero  valued  outside  a  radius  of  five  nodes  from  the 
origin.  The  domain  is  comprised  of  73  equi-spaced  nodes  in  each  direction;  non-reflecting 
boundary  conditions  were  applied  on  all  six  sides  of  the  domain.  Three  points  in  space  will 
be  examined:  the  origin,  an  arbitrary  point  inside  the  initial  perturbation  (r  <  5 dx),  and 
an  arbitrary  point  outside  the  initial  perturbation  (r  >  5 dx).  Referring  to  the  notation  of 
Equations  (55)  -  (58), 


<p(x,y,z) 


1  —  -  if  r  <  a, 
0  if  r  >  a 


(71) 


^(x,y,z)  =  0 


(72) 


where  r  is  the  cartesian  radius,  \J  x2  +  y2  +  z2.  The  analytic  solution  of  the  integrals  in 
Equation  (58)  is  not  easily  calculated,  but  application  of  a  composite  Simpson’s  Rule  over 
intervals  of  ^  will  yield  a  suitable  comparison. 

Figures  48  -  50  show  the  analytical  solutions  at  the  respective  points  plotted  as  a 
function  of  time  directly  over  their  CA  counterparts  plotted  versus  the  number  of  iterations 
through  which  the  rule  has  been  applied.  These  plots  illustrate  a  key  challenge  to  CA  as 
identified  by  Hosoglu  [46]:  matching  the  discrete  iterations  of  a  cellular  automaton  to  the 
continuous  time  domain,  or  calculating  an  appropriate  time  scaling  factor  (TSF).  In  one 
dimension,  dt  =  —  works  well,  but  for  three  dimensions,  as  Figure  51  shows,  there  is  a 
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phase  difference  resulting  from  a  time  scale  mismatch.  Consider  a  true  point  source  located 
at  the  origin  of  an  otherwise  zero-valued  CA  domain.  Following  the  current  CA  rule,  the 
earliest  a  node  at  (dx,  dy,  dz )  can  reach  a  value  other  than  zero  is  after  the  third  iteration; 
thus  3 dt  =  -c\J dx2  +  dy 2  +  dz2,  or,  for  an  equi-axed  mesh,  dt  =  ^===.  The  points  of 
interest  for  this  exercise  were  chosen  arbitrarily,  but  in  such  a  way  that  both  the  inside  and 
outside  points  are  displaced  from  the  origin  in  all  three  directions.  Comparisons  between 
the  analytical  and  CA  solutions  plotted  with  this  time  equivalency  are  shown  in  Figures 
52-54. 


analytic  solution  at  (0,0,0) 

1 - - - . - - - 


0  12  time  3  4  5 


Figure  48:  Three-dimensional  wave  model  at  domain  origin:  time  vs.  iterations 
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analytic  solution  at  (-dx,  -3dx,  -2dx)  [r=3.75dx] 
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Figure  49:  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturbation:  time 
vs.  iterations 
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Figure  50:  Three-dimensional  wave  model  at  a  point  outside  the  initial  perturbation:  time 
vs.  iterations 
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time  scale  mis-match 


Figure  5 1 :  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturbation:  analytic 
solution  vs  time,  CA  solution  vs  — 

’  C 


origin 


Figure  52:  Three-dimensional  wave  model  at  a  domain  origin 
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inside 


Figure  53:  Three-dimensional  wave  model  at  a  point  inside  the  initial  perturbation 


outside 


Figure  54:  Three-dimensional  wave  model  at  a  point  outside  the  initial  perturbation 
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D.  COUPLING  OF  FINITE  ELEMENT  AND  CELLULAR  AUTOMATA  MOD¬ 
ELS  OF  THE  WAVE  EQUATION 

The  CA  model  of  the  wave  equation  for  velocity  potential  is  simple  to  implement 
and  can  easily  be  adapted  to  a  variety  of  non-trivial  boundary  conditions,  but  converting 
velocity  potential  back  into  useful  quantities  like  pressure  and  velocity  as  time-  and  spatial- 
derivatives  is  hampered  by  the  alternating  update  nature  of  the  model.  The  finite  element 
model,  on  the  other  hand  updates  the  value  of  every  point  every  time-step,  but  can  be  com¬ 
putationally  expensive  and  non-trivial  boundary  conditions  can  be  difficult  to  implement. 
A  combination  of  the  two  methods  would  seem  to  resolve  the  short-comings  of  each  and 
enable  a  faithful  model  of  the  fluid  portion  of  our  fluid-structure  interaction. 

The  general  idea  is  to  have  several  layers  of  finite  elements  in  contact  with  the 
structure  and  then  to  have  that  fluid  volume  surrounded  with  a  CA  domain  upon  which  the 
non-reflecting  and  free  boundary  conditions  can  be  imposed  as  shown  in  Figure  55.  The 
two  fluid  domains  will  overlap  such  that  the  outer  layer  of  finite  element  nodes  will  be  pro¬ 
cessed  as  interior  CA  nodes  whose  CA-calculated  values  become  FE-specified  boundary 
values.  The  next  set  of  FE  nodes  inside  the  domain  are  calculated  by  the  FE  machinery  and 
then  passed  to  the  CA  domain  to  serve  as  neighbors  for  application  of  the  CA  rule  to  the 
outer  set.  These  node  sets  are  illustrated  in  Figure  56. 

The  pre- validation  of  this  scheme  was  to  establish  comparable  domains  of  each 
model,  impose  a  specified  velocity  potential  field  on  one  face  of  the  domains  (the  top), 
and  specify  a  fixed,  zero-valued  boundary  condition  on  the  other  five  faces.  The  specified 
input  is  a  radially  scaled  sinusoid-it  achieves  its  maximum  value  at  the  center  of  the  face 
over  which  it  is  applied  and  is  zero  beyond  a  radius  of  one-half  the  width  of  the  region. 
The  resulting  value  of  the  velocity  potential  at  the  respective  domain  centers  compares 
favorably  as  shown  in  Figure  57. 

The  first  full  validation  of  the  coupling  scheme  was  to  model  a  joint  domain  with 
an  imposed  velocity  potential  on  the  top  of  the  FE  portion  of  the  domain  which  in  turn 
rests  atop  the  CA  portion.  The  four  sides  of  both  domains  have  fixed  zero-valued  boundary 
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finite  element  domain  within  cellular  automata  domain 


y  axis 


x  axis 


Figure  55:  Finite  Element  fluid  domain  surrounded  by  Cellular  Automata  fluid  domain 

conditions  and  the  bottom  of  the  CA  portion  is  non-reflecting.  The  value  of  the  velocity  po¬ 
tential  on  either  side  of  the  interface  is  shown  in  Figure  58.  Next,  the  specified  function  on 
the  top  of  the  FE  domain  (the  same  radially  scaled  sinusoid)  was  treated  as  a  velocity  rather 
than  velocity  potential.  Again,  the  values  of  0  are  compared  near  the  interface  between  the 
two  models  and  shown  in  Figure  59.  Next,  the  domain  and  interfaces  are  expanded  such 
that  the  CA  domain  surrounds  the  FE  domain  on  five  sides  and  the  non-reflecting  boundary 
condition  is  applied  to  all  six  sides  of  the  outer  domain  with  the  exception  of  the  top  of  the 
FE  domain-velocity  is  specified  there.  Those  results  are  shown  in  Figure  60. 

The  time-integration  of  Equation  (65)  was  performed  using  a  New  mark-/)  algorithm 
with  a  zero  valued  [C]  matrix. [47]  The  cases  that  involved  fixed,  zero  boundary  conditions 
(Figures  57  -  59)  display  noticeable  high-frequency  noise  that  is  suspected  to  be  caused 
by  the  sudden  change  in  <f>  value  at  the  boundaries.  Despite  a  smooth  input  function,  the 
beginning  of  a  similar  phenomenon  is  noticeable  in  the  case  with  completely  non-reflecting 
boundary  conditions  as  well  (Figure  60).  The  particular  Ncwmark-/)  scheme  used  (7  =  |, 
/3  =  |)  is  unconditionally  stable  and  was  chosen  to  ameliorate  any  potential  difficulty 
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O  FE  calculation  passed  to  CA 
+  CA  calculation  passed  to  FE 


Figure  56:  Node  sets  for  exchange  of  data  between  fluid  domains 


arising  from  dt  being  a  fixed  function  of  dx,  however  a  transition  to  an  o-mcthod  (or 
Hilber-Hughes-Taylor  (HHT))  may  be  necessary  to  dampen  this  noise.  See  Appendix  A 
for  details  on  these  methods. 
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velocity  potential  at  center  of  domain 


Figure  57:  Comparison  of  velocity  potential  propagation  between  finite  element  and  cellu¬ 
lar  automata  models  with  common  Dirichlet  boundary  conditions 


velocity  potential  near  (0,0, z),  <|>(x,y,^o  ,t)  specified 


Figure  58:  Velocity  potential  propagation  between  finite  element  and  cellular  automata 
domains,  velocity  potential  (©(top))  specified 
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velocity  potential  near  (0,0, z.),  v(x,y,^op,t)  specified 


Figure  59:  Velocity  potential  propagation  between  finite  element  and  cellular  automata 
domains,  velocity  (w(top))  specified 


velocity  potential  near  (0,0, z.),  v(x,y,^op,t)  specified 


Figure  60:  Velocity  potential  propagation  between  finite  element  and  cellular  automata 
domains,  FE  inside  CA,  velocity  (u(top))  specified,  non-reflecting  boundary  conditions 
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1.  Comparison  with  Homogeneous  Fluid  Domain 

Figure  61  demonstrates  close  agreement  between  the  response  of  a  composite  do¬ 
main  and  a  homogeneous  fluid  domain  subject  to  the  same  input  and  boundary  conditions. 


velocity  potential  at  mid-domain 


Figure  61:  Comparison  of  velocity  potential  at  mid-domain  resulting  from  specified  value 
on  one  face:  FE  +  CA  domain  vs.  homogeneous  CA  domain 
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V.  RESULTS 


A.  ACOUSTIC  FIELD  FLUID-STRUCTURE  INTERACTION 

This  chapter  will  combine  the  structural  and  fluid  models  discussed  so  far  and  ex¬ 
amine  their  interactions.  Finally,  comparison  with  experimental  results  will  be  presented. 

1.  Information  Exchange 

As  previously  discussed,  the  fluid  model  used  in  this  work  is  that  of  the  velocity 
potential.  At  each  time  step  the  transverse  velocity  of  the  structure  is  transmitted  into  the 
fluid  where  it  is  converted  to  velocity  potential  and  propagated  through  the  fluid  domain  ac¬ 
cording  to  the  wave  equation.  The  resulting  fluid  pressure,  a  function  of  the  time  derivative 
of  the  velocity  potential  is  then  applied  to  the  structure  via  its  load  vector  for  calculation 
of  displacement  and  velocity  during  the  next  time  step.  Unconditionally  stable  time  inte¬ 
grators  across  the  various  domains  allows  the  selection  of  time  step  size  based  on  the  CA 
time-scaling  factor. 

For  simplicity  the  meshes  of  the  structure,  the  FE  fluid,  and  the  CA  fluid  are  mutu¬ 
ally  conforming. 

2.  Homogeneous  Isotropic  Single-Layer  Structure 

The  algorithm  described  above  was  validated  by  modeling  a  clamped  foot  square 
plate  one-quarter  inch  thick  homogeneously  comprised  of  an  isotropic  material.  Both  CG 
and  DG  models  were  generated.  Fresh  water  material  properties  were  used  for  the  fluid 
model.  Figure  62  shows  the  transverse  displacement  of  the  center  of  the  plate  subject  to 
a  constant  concentrated  force  applied  at  its  center.  The  dry  case  is  the  resulting  oscil¬ 
lation  about  its  predicted  static  deflection;  the  wet  case  shows  that  both  magnitude  and 
frequency  of  the  oscillation  have  been  altered  as  a  result  of  the  FSI.  Kwon  [48]  discussed 
this  phenomena  and  observed  the  effects  of  the  elastic  modulus  and  density  of  the  struc¬ 
ture.  Specifically,  he  noted  that  a  structure  with  a  density  close  to  that  of  the  fluid  would 
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be  more  affected  by  the  interaction.  This  observation  is  borne  out  in  Figures  63  and  64  for 
which  the  respective  material  properties  were  changed  from  the  baseline  case.  The  baseline 
model  has  a  structural  density  2.7  times  that  of  the  fluid,  a  frequency  ratio  (dry /wet)  of  1 .84 
and  an  amplitude  ratio  (first  peak  -  first  trough,  dry/wet)  of  2.49.  The  double  modulus  case 
(Figure  63)  shows  a  slightly  higher  frequency  oscillation  about  a  lesser  static  deflection  in 
the  dry  case;  the  dry/wet  frequency  ratio  is  1 .86  and  amplitude  ratio  is  2.15:  a  lesser  relative 
change.  The  double  density  case  (Figure  64)  shows  the  expected  lower  frequency  oscilla¬ 
tion  in  the  dry  case,  but  also  a  dry /wet  frequency  ratio  of  1.49.  The  dry /wet  amplitude  ratio 
for  the  double  density  case  is  2.00. 


Figure  62:  Displacement  of  clamped  plate  with  and  without  fluid- structure  interaction 
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Figure  63:  Displacement  of  clamped  plate  of  double  modulus  with  and  without  fluid 
structure  interaction 


Figure  64:  Displacement  of  clamped  plate  of  double  density  with  and  without  fluid 
structure  interaction 
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3.  Two-layer  Plates 

The  clamped  plates  modeled  in  this  section  are  each  0.3048m  x  0.3048m  x  3.5mm 
and  comprised  of  two  thickness  layers  and  sixteen  elements  in  each  planar  direction.  The 
material  is  an  isotropic  approximation  of  E-glass.  Initial  conditions  were  zero  displacement 
and  velocity;  a  constant  concentrated  force  of  1000N  applied  at  the  center  of  the  plate  at 
the  first  time  step.  The  “damaged”  plates  had  a  four  element  by  four  element  debonding 
patch  inserted  between  the  two  layers  at  the  center  of  the  plate.  This  debond  was  of  the 
partial  disconnection  method  described  in  Chapter  III.  Figures  65  -  67  display  the  time 
histories  of  displacement,  velocity,  and  normal  strain  in  the  plane  on  the  bottom  of  each 
plate.  The  stress  profiles  in  Chapter  III  showed  that  maximum  values  were  at  the  edges  of 
the  damage  zone;  similar  phenomena  are  expected  in  strain  values,  but  as  Figure  68  shows, 
the  strains  at  the  center  of  the  plate  for  both  dry  and  wet  cases  is  of  greater  magnitude  than 
those  at  the  +y  edge  of  the  damage  zone  (the  maximum  numerically  of  the  four  edges). 
This  further  suggests  that  an  interface  layer  is  necessary  to  properly  model  debonding  in 
laminated  composites  and  other  layered  structures.  The  density  of  E-glass  is  twice  that  of 
water,  close  enough  to  observe  an  added  mass  effect  in  the  case  of  FSI.  That  is,  a  structure 
with  density  comparable  to  that  of  fluid  reacts  not  only  to  the  pressure  effect  of  the  fluid- 
dampening  the  amplitude  of  its  vibrations,  but  also  at  a  lower  frequency  as  though  it  was  a 
more  massive  structure. 
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damaged  two  layer  E-glass  plate  displacement 


Figure  65:  Displacement  of  damaged  clamped  two  layer  E-glass  plate  with  and  without 
fluid-structure  interaction 


damaged  two  layer  E-glass  plate  velocity 


Figure  66:  Velocity  of  damaged  clamped  plate  two  layer  E-glass  plate  with  and  without 
fluid-structure  interaction 
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Figure  67:  Strain  at  center  of  clamped  two  layer  E-glass  plate  with  and  without  fluid- 
structure  interaction 


Figure  68:  Strain  at  center  and  edge  of  damage  zone  of  clamped  two  layer  E-glass  plate 
with  and  without  fluid- structure  interaction 
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4.  Three-layer  Plates 

The  two-layer  model  of  the  previous  sub-section  did  a  poor  job  of  reflecting  debond¬ 
ing  damage  within  a  laminated  plate.  A  three-layer  model  with  a  thin  interface  layer  with 
properties  approximating  a  common  adhesive  inserted  between  two  layers  of  E-glass  was 
subjected  to  the  same  loading  and  boundary  conditions  (zero  initial  displacement  and  ve¬ 
locity,  1000  N  concentrated  force  at  center,  clamped  edges);  responses  were  calculated  for 
five-hundred  time  steps.  Figures  69  and  70  show  that  the  displacement  and  velocity  re¬ 
sponses  of  the  center  of  the  plate  do  not  reflect  the  presence  or  absence  of  a  debonding 
zone,  but  do  demonstrate  FSI  mass  effects  in  a  fashion  similar  to  that  of  the  two-layer 
model.  Figure  71  shows  that  the  strain  calculated  in  the  E-glass  elements  reflects  presence 
or  absence  of  damage  only  mildly.  Figure  72,  on  the  other  hand,  shows  clearly  that  the 
interface  layer  is  profoundly  affected  by  the  presence  of  a  damage  zone.  To  be  clear,  the 
strain  values  at  the  centers  of  the  interface  layers  of  the  damaged  plates  are  not  identically 
zero,  but  they  are  four  orders  of  magnitude  lower  than  their  undamaged  counterparts. 


three  layer  E-glass  plate  displacement 


Figure  69:  Displacement  of  clamped  three  layer  E-glass  plate  with  and  without  fluid- 
structure  interaction  and  with  and  without  damage 
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three  layer  E-glass  plate  velocity 


Figure  70:  Velocity  of  clamped  three  layer  E-glass  plate  with  and  without  fluid-structure 
interaction  and  with  and  without  damage 


Figure  71:  Strain  of  clamped  three  layer  E-glass  plate  with  and  without  fluid-structure 
interaction  and  with  and  without  damage  at  center  of  lower  E-glass  layer 
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Figure  72:  Strain  of  clamped  three  layer  E-glass  plate  with  and  without  fluid-structure 
interaction  and  with  and  without  damage  at  center  of  interface  layer 
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Figure  73:  Strains  in  dry  clamped  three  layer  E-glass  plate  with  and  without  damage 
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Figure  74:  Strains  in  wet  clamped  three  layer  E-glass  plate  with  and  without  damage 
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B.  EXPERIMENTAL  VALIDATION 


Recent  experimental  work  by  Kwon  and  his  students  [49],  [50],  [51],  [29]  examined 
the  response  of  composite  plates  to  low  velocity  impact  with  and  without  fluid- structure  in¬ 
teraction^).  In  general,  they  found  that  for  a  given  impact  weight  dropped  from  the  same 
height,  structures  with  lower  density  relative  to  the  fluid  in  question  experienced  higher 
resultant  forces  and  consequently  greater  damage  than  the  same  material  in  dry  conditions. 
Additionally,  they  noted  that  the  initial  observable  damage  mode  was  delamination  occur¬ 
ring  on  the  face  of  the  plate  opposite  the  impact  site. 
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Figure  75:  Schematic  of  VARTM  for  plate  manufacture.  From  [49] 


In  his  work,  Conner  [29]  used  a  Vacuum  Assisted  Resin  Transfer  Molding  technique, 
shown  in  Figure  75,  to  construct  a  series  of  12in  by  12in  composite  plates  comprised  of 
sixteen  layers  of  E-glass  (approximately  3.5mm  thick  in  toto)  and  subjected  them  to  low- 
velocity  impact  forces  that  resulted  from  dropping  a  10.8kg  weight  from  various  heights  to 
the  center  of  the  plate(s)  using  the  assembly  shown  in  Figures  76  and  77.  The  plates  were 
instrumented  with  strain  rosettes  at  four  set  positions  and  oriented  such  that  one  channel 
returned  ex  directly;  ey  was  calculated  as  a  function  of  all  three  channels  and  the  included 
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Figure  76:  Drop  Weight  Rig  used  in  Impact  Testing.  From  [48] 

angles  as  can  be  found  in  many  solid  mechanics  textbooks  including  [52].  Data  was  sam¬ 
pled  at  a  frequency  of  10,000  Hz  (dt  =  10-4s). 

Numerical  comparison  with  this  experimental  data  was  conducted  using  a  DG  struc¬ 
tural  model  comprised  of  a  single  layer  of  plate  elements  with  a  discretization  of  twelve 
elements  in  each  planar  direction.  The  overall  structure  has  a  length  to  thickness  ratio  of 
87:1  and  each  element  has  a  length  to  thickness  ratio  of  7.3:1.  In  this  model  the  material 
properties  used  are  those  of  E-glass,  but  treated  as  an  isotropic  material-the  Young’s  mod¬ 
ulus  along  its  fiber  direction  was  taken  in  all  three  directions.  The  mass  matrix  is  lumped 
and  therefore,  diagonal.  The  dry  response  was  calculated  using  the  a-method  time  inte¬ 
gration  of  the  equation  of  motion  for  the  plate.  The  wet,  or  FSI,  response  was  calculated 
according  to  the  acoustic  field  FSI  described  in  the  last  chapter  with  time  step  size  for  the 
entire  model  equal  to  the  TSF  for  the  CA  portion  of  the  fluid  model. 

The  force  inputs  to  the  numerical  plates  were  smoothed  versions  of  the  experimen¬ 
tally  measured  force  data  in  the  time  region  of  interest-the  main  impact.  The  smoothing 
was  conducted  by  sampling  the  raw  data  every  five  time  steps,  generating  an  interpolation 
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Figure  77:  Drop  Weight  Rig  as  used  for  Impact  Testing  with  FSI.  From  [48] 

function  (using  MATLAB’s  interpl  function  with  the  spline  option)  and  evaluating 
the  interpolation  function  at  every  required  time  step  for  the  model.  For  the  dry  cases,  the 
same  time  step  size  as  the  experimental  was  used.  Figure  78  shows  the  raw  experimental 
force  data  and  the  smoothed  version  as  generated  for  the  dry  case;  Figure  82  shows  the 
force  data  used  for  the  FSI  case. 

The  calculated  time  history  of  the  displacement  field  was  used  to  calculate  a  time 
history  of  the  strain  vector  at  nodal  points  throughout  the  domain  of  the  plate.  Those  nodes 
closest  to  the  positions  of  the  strain  gages  in  the  experimental  work  were  examined  relative 
to  the  recorded  data  and  are  plotted  in  Figures  79-81  for  the  dry  plate  and  Figures  83  -  85 
for  the  FSI  case.  All  show  good  qualitative  agreement  between  experimental  and  numerical 
data.  Differences  can  be  attributed  to  the  smoothing  of  the  input  force,  ignoring  impact 
effects,  the  isotropic  treatment  of  an  orthotropic  material,  the  approximation  of  structural 
thickness,  approximate  positioning  of  strain  gages,  and  mis-alignment  of  strain  gages.  The 
data  correlating  with  gage  1  appears  to  show  better  overall  agreement  than  the  other  two, 
most  likely  because  that  gage  was  located  approximately  equi-distant  from  both  the  point 
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of  impact  and  the  clamped  boundaries  of  the  plate.  The  other  gages  were  closer  to,  and 
therefore,  more  exposed  to  the  effects  of  the  physical  boundaries  of  the  plate. 


force  input 


Figure  78:  Raw  and  smoothed  experimental  force  data  for  dry  plate 
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Figure  79:  Measured  versus  calculated  strain,  dry  plate,  gage  1 


ex  gage  2  £y  gage  2 


Figure  80:  Measured  versus  calculated  strain,  dry  plate,  gage  2 
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Figure  81:  Measured  versus  calculated  strain,  dry  plate,  gage  3 


force  input 


Figure  82:  Raw  and  smoothed  experimental  force  data  for  wet  plate 
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Figure  83:  Measured  versus  calculated  strain,  wet  plate,  gage  1 


ex  gage  2  ey  9a9e  2 


Figure  84:  Measured  versus  calculated  strain,  wet  plate,  gage  2 
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Figure  85:  Measured  versus  calculated  strain,  wet  plate,  gage  3 
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VI.  CONCLUSION 


A.  SUMMARY  OF  FINDINGS 

The  goal  of  this  work  was  to  develop  computational  techniques  to  accurately  model 
and  simulate  dynamic  responses  and  failures  of  composite  structures  in  an  acoustic  field. 
After  implementing  a  nodal  three-dimensional  element  to  verify  basic  computational  method¬ 
ology,  a  displacement-only  plate  finite  element  was  formulated  and  implemented  using 
Discontinuous  Galerkin  (DG)  methodology.  Such  a  displacement-only  element  allows  con¬ 
struction  of  multi-layered  structures  like  sandwich  plates  and  other  laminated  composites 
in  a  manner  similar  to  full  three-dimensional  solid  finite  elements.  Results  generated  from 
this  formulation  compare  favorably  with  theoretical  predictions  as  well  as  existing  CG  nu¬ 
merical  models  for  both  static  and  dynamic  responses  for  both  simple  and  multi-layered 
plate  structures. 

Application  of  the  new  element  to  the  analysis  of  failure  initiation  and  propagation 
in  sandwich  composite  structures  shows  great  promise.  Static  qualitative  stress  profiles  are 
similar  to  those  found  using  CG  techniques,  but  the  elemental  rather  than  nodal  connectivity 
used  in  DG  formulations  suggests  a  simple  means  of  modeling  debonding  between  material 
layers  by  disconnecting  their  respective  elements  in  the  global  stiffness  matrix.  Complete 
disconnection  of  neighboring  elements  in  an  imposed  debonding  zone  was  shown  to  be 
incorrect  because  it  allowed  the  core  layer  to  deflect  not  only  through  the  disconnected 
resin  layer  but  also  through  the  still-present  skin  layer.  Partial  disconnection-removing 
connectivity  between  opposing  pairs  of  dofs  in  the  planar  directions  but  retaining  weak 
connectivity  for  transverse  pairs  of  dofs  resulted  in  a  stress  profile  that  makes  good  quali¬ 
tative  sense.  Maximum  stress  values  in  both  skin  and  core  layers  decreased  in  magnitude 
and  moved  from  the  center  of  the  plate  to  the  edges  of  the  debonding  zone,  behaving  like  a 
stress  concentration.  The  simplicity  of  this  partial  disconnection  method  can  be  a  tremen- 
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dous  computational  savings  when  modeling  the  progression  of  damage  without  need  for 
re-meshing  or  recalculation  of  the  global  stiffness  matrix. 

Examination  of  FSI  for  the  impact  problem  does  not  require  a  full  fluid-flow  model- 
the  propagation  of  velocity  potential  according  to  the  wave  equation  in  the  acoustic  domain 
is  sufficient  for  these  purposes;  as  such,  an  extension  of  Cellular  Automata  (CA)  from  two 
to  three-dimensions  in  modeling  the  acoustic  field  was  demonstrated  and  validated.  CA 
was  chosen  for  this  application  not  only  because  of  the  simplicity  of  its  update  rule,  but 
also  because  of  its  flexibility  in  the  implementation  of  non-reflecting  boundary  conditions. 
The  alternating  update  nature  of  CA  makes  calculation  of  both  spatial  and  time  derivatives- 
required  to  convert  the  velocity  of  the  structure  into  velocity  potential  in  the  fluid  domain 
and  to  convert  velocity  potential  into  a  pressure  field-difficult.  Insertion  of  a  small  finite  el¬ 
ement  interface  zone  between  the  structure  and  the  CA  fluid  domain  resolved  this  difficulty 
for  a  relatively  low  computational  cost.  The  combination  of  a  small  FE  acoustic  domain 
with  an  enveloping  CA  domain  proved  to  be  an  efficient  way  to  implement  non-reflecting 
boundary  conditions. 

Finally,  the  combined  model  of  a  DG  structure  interacting  with  a  FE+CA  fluid  do¬ 
main  was  shown  to  have  good  agreement  between  calculated  and  experimentally  measured 
strain  values  for  plates  subject  to  low-velocity  impact  in  the  pre-damage  regime.  In  partic¬ 
ular,  the  added  mass  effect  on  structures  with  low  density  relative  to  the  fluid  medium  was 
apparent  in  both  simulation  and  experimental  comparisons. 

The  methods  developed  and  examined  in  this  study:  the  displacement-only  DG 
plate  finite  element,  the  partial  disconnection  failure  model,  and  the  hybrid  FE+CA  acoustic 
field  model,  show  great  promise  for  flexible  and  accurate  modeling  of  debonding  initiation 
and  propagation  in  sandwich  and  laminate  composite  structures  subject  to  FSI. 

B.  FUTURE  WORK 

This  work  should  be  viewed  as  a  starting  point  for  further  investigation  into  the 
utility  of  DG  methods  in  composite  failure  modeling. 
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First  and  foremost,  the  current  formulation  must  be  re-implemented  to  be  a  true 
element-wise  computation  in  order  to  reap  the  benefits  of  DG  not  just  pay  the  costs  of  solv¬ 
ing  for  a  greater  number  of  dof.  This  type  of  update  will  enable  more  efficient  and  flexible 
modeling  of  larger  problems,  including  more  refined  meshes,  more  complex  geometry,  and 
approaches  that  address  all  levels  of  multi-scale  analysis  of  composite  materials.  Such  a 
re-implementation  should  also  enable  the  coupling  of  the  current  DG  formulation  with  CG 
codes.  This  should  be  readily  achievable  as  the  current  formulation  is  derived  from  one 
with  such  coupling  as  a  specific  goal.  This  sort  of  coupling  can  be  used  as  an  alternative  to 
refining  the  mesh  in  the  areas  of  interest  like  existing  or  expected  damage  by  replacing  the 
refined  mesh  with  DG  elements  in  order  to  better  examine  the  physical  phenomena. 

A  more  computationally  efficient  implementation  should  also  include  and  enable 
progressive  failure  modeling  both  through  inclusion  of  traction-separation  type  models  for 
the  post-damage  regime,  but  also  through  propagation  of  damage  beyond  its  initiation  site. 
Addition  of  a  full  impact-impulse  model  for  force  input  should  enable  even  more  faithful 
modeling  and  closer  comparison  with  experimental  results. 

The  closer  a  computational  model  approaches  observed  physical  phenomena,  the 
more  useful  and  trustworthy  its  results  in  evaluating  more  complex  geometries  and  operat¬ 
ing  environments. 


97 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


98 


APPENDIX  A.  TIME  INTEGRATION  ALGORITHMS 


This  appendix  contains  explication  of  the  algorithms  used  to  solve  the  matrix- vector 
equation  of  motion  in  this  work.  Implicit  methods  are  favored  for  their  stability  independent 
of  size  of  time  step-a  concern  when  matching  various  domains.  All  methods  are  trying  to 
solve 

[M]{u}  +  [C]{u}  +  [A-]{tl}  =  {/}  (73) 

for  {u}.  Solutions  for  the  two  time  derivatives  are  used  as  needed  to  update  {»};  in  the 
finite  element  fluid  domain  {«}  is  also  used  for  a  pressure  calculation. 

A.  NEWMARK-/3  METHOD 

The  below  algorithm,  taken  from  [47]  was  initially  implemented  to  serve  as  time 
integrator  of  Equation  (73)  by  using  the  weighted  averages 

{u}n+ 1  =  {u}n  +  [(1  -  l){u}n  +  l{u}n+l]  '  dt  (74) 

(fa2 

{u}n+ 1  =  {u}n  +  {u}n  ■  dt  +  [(1  -  2 fd){u}n  +  2/3{u}n+1]  ■  —  (75) 

substituting  and  rearranging  terms  results  in 

[M  +  ^dtC  +  pdt2K]{u}n+1  =  {f}n+ 1 

dt2 

-  [(1  -  7 )dtC  +  (1  -  2 (3)—K]{u}n  -  [C  +  dtK]{u}n  -  K{u}n  (76) 

which  can  be  solved  for  u  and  then  u  and,  in  turn,  u  at  each  time  step.  The  method  is 
unconditionally  stable  for  2/3  >  7  >  |.  Parameter  choices  are  7  =  \  and  /3  —  \  correspond 
to  Newmark’s  Constant  Average  Acceleration  Method. 

Dirichlet  and  Neumann  boundary  conditions  can  be  imposed  nodally  by  solving 
for  the  current  acceleration  on  the  boundary  through  Equations  (74)  and  (75),  zeroing  the 
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corresponding  rows  of  the  compound  left-hand-side  matrix,  setting  the  diagonal  elements 
of  those  rows  to  1  and  substituting  the  boundary  accelerations  into  the  right-hand-side 
vector. 

B.  q-METHOD 

Hilber,  Hughes,  and  Taylor  improved  upon  Newmark-/?  with  their  introduction  of 
the  a-method  [53]-[54].  This  method  is  designed  to  dissipate  high  frequency  noise  without 
degrading  the  order  of  solution  accuracy.  The  update  rules  for  {  «}  and  {  «}  are  the  same  as 
in  (75)  and  (74),  but  now  the  equation  of  motion  is  also  a  weighted  average: 

[M]{u}n+ 1  +  (1  +  a)[C]{ii}n+l  -  Oi[C]{u}n 

+  (1  +  a)[K]{u}n+ 1  -  a[K}{u}n  =  {f{tn+ i+a)}  (77) 


which  is  rearranged  to 

[M  +  (1  +  a)dt{jC  +  /3dtK)\{u}n+1  =  (1  +  a)fn+1  -  afn 

-(l+a)dt[(\-1)C+d^(l-2p)K){u}n 

~[C  +  dt(  1  +  a)K]{u}n  -  K{u}n  (78) 

and  then  solved  for  {«}n+1  which  is  then  used  to  update  {«}n+1,  and  (w}n+i.  Dirichlet 
boundary  conditions  are  applied  nodally  via  rearrangement  of  Equation  (75)  to  solve  for 
prescribed  values  of  the  right  hand  side  of  Equation  (78).  This  method  is  unconditionally 
stable  when  a  £  [— |,  0],  7  =  (1  —  2a)/2,  and  —  (1  —  a)2/ 4.  When  a  =  0  this  method 
reduces  to  Newmark’s  Constant  Average  Acceleration  Method. 

C.  TIME  DISCONTINUOUS  GALERKIN  METHOD 

Another  time  integrator  considered  but  not  fully  implemented  in  this  work  is  the 
Time  Discontinuous  Galerkin  (TDG)  method  presented  by  Chien,  Yang,  and  Tang  [55]. 
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They  present  time  as  yet  another  domain  that  can  be  discretized  by  finite  elements,  in  par¬ 
ticular  as  discontinuous  finite  elements,  as  shown  in  Figure  86.  For  the  undamped  equation 
of  motion  ([C]=0  in  Equation  (73)),  they  generate  the  following  matrix  equation  to  solve 
for  the  displacement  and  velocity  at  each  end  of  a  particular  time  interval  (or  element): 
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where  K  and  M  are  the  usual  stiffness  and  mass  matrices  and 


F\  = 


0i(f)Fdf 


(80) 


F2=  [  fa(t)Fdt  (81) 

Jin 

where  F  is  the  usual  applied  load  vector  and  0\ (t)  and  02 (t)  are  the  time  shape  functions 
shown  in  Figure  86.  Simplifying  and  recasting  some  terms  results  in 
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which  can  be  solved  for  the  velocity  terms  from  which  displacements  may  be  calculated 
directly.  They  also  show  that  this  method  is  also  unconditionally  stable,  making  it  an 
intriguing  avenue  for  future  work  in  structural  dynamics. 


Figure  86:  Temporal  elements  for  TDG  method.  From  [55] 
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APPENDIX  B.  IMPLEMENTATION  DETAILS 


The  discussion  of  the  fluid  and  structural  and  fluid  models  in  the  main  body  of  this 
thesis  is  (moderately)  general  and  symbolic.  All  implementation  was  in  MATLAB,  and  run 
on  a  MacBook  Air  (mid-2009),  MacBook  Pro  (mid-201 1),  or  on  the  Hamming  cluster. 

A.  MESHING 

The  various  domains  were  meshed  for  human  convenience  rather  than  any  matrix 
bandwidth  considerations.  In  general,  node  and  element  numbering  proceeds  from  (xmin, 
Dmin ,  Zmin )  along  the  x-axis,  then  increment  up  the  y-axis,  and  then  up  the  z- axis.  While  the 
code  does  execute  full  and  proper  Jacobian  calculations  for  arbitrarily  oriented  elements,  in 
general  the  x-axis  corresponds  to  the  canonical  r-axis;  the  y-axis  corresponds  to  the  canon¬ 
ical  s-axis,  and  the  x-axis,  which  also  generally  corresponds  to  thickness,  corresponds  to 
the  canonical  t-axis. 

The  various  meshes  employed  are  conforming  for  human  convenience.  The  fluid 
domain  is  entirely  equi-axed.  The  FE  portion  is  constructed  using  the  coordinates  of  the 
bottom  of  the  plate  elements  as  a  foundation  and  extending  down  a  specified  number  of 
layers  in  the  — z  direction.  The  CA  portion  is  constructed  by  specifying  a  factor  by  which 
plate  length  is  multiplied-the  cube  of  this  value  is  the  fluid  volume.  The  CA  nodes  that  are 
wholly  inside  the  FE  fluid  domain  and  not  needed  are  simply  removed  from  the  index  sets 
and  ignored. 

B.  NUMERICAL  INTEGRATION 

Exact  integration  is  performed  using  Gauss-Lobatto  quadrature.  Under-integration 
of  the  shear  terms  of  plate  elements  is  performed  using  Gauss-Legendre  quadrature.  Be¬ 
cause  the  integration  points  for  Gauss-Lobatto  quadrature  are  still  nodal  (interpolation) 
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points,  terms  in  Ks  that  are  products  of  different  cardinal  interpolation  functions  will  be 
uniformly  zero,  resulting  in  a  too  sparse  shear  stiffness  matrix. 

A  numerical  experiment  comparing  the  static  deflections  calculated  using  the  two 
different  quadrature  rules  for  calculating  an  under- integrated  K s  with  both  theoretical  static 
deflection  and  that  resulting  from  an  exactly  integrated  Ks  was  conducted.  For  a  clamped 
plate  of  dimensions  0.3048m  x  0.3048m  x  0.00635m  (12in  x  12in  x  l/4in),  elastic  modulus 
70  GPa,  Poisson’s  ratio  0.3  subjected  to  a  concentrated  load  of  1000N,  predicted  static 
deflection  of  the  center  of  the  plate  is  0.31697mm  [33].  Calculated  deflection  for  Gauss- 
Lobatto  quadrature  was  0.0945mm;  for  Gauss-Legendre  quadrature:  0.3086mm;  for  exact 
integration  of  Ks:  0.0325mm.  Clearly,  Gauss-Lengendre  quadrature  is  a  better  choice  for 
under-integration  of  the  shear  stiffness  matrix  in  this  application. 

C.  APPLICATION  OF  BOUNDARY  CONDITIONS  AND  EXTERNAL  LOADS 

Boundary  conditions  and  external  loads  are  applied  nodally.  For  boundary  condi¬ 
tions,  rows  of  the  mass  and  stiffness  matrices  corresponding  to  constrained  dofs  are  zeroed; 
their  diagonal  elements  are  set  to  1,  and  the  corresponding  element  in  the  right-hand-side 
vector  is  set  to  the  constrained  value.  External  loads  must  first  be  converted  to  units  of 
force  and  then  distributed  to  the  appropriate  elements  of  the  right-hand-side  vector.  For 
the  pressure  exerted  by  the  fluid  domain  on  the  wet  side  of  our  notional  plate,  the  nodal 
pressure  vector  is  multiplied  by  a  two-dimensional  interpolation  (mass)  matrix  to  convert  it 
to  a  force  vector.  For  application  to  the  DG  structure,  that  resulting  vector  is  decomposed 
to  reflect  the  number  of  dof  found  at  each  geometric  position;  in  this  way,  the  force  applied 
at  a  point  ’’shared”  by  four  separate  discontinuous  elements  will  be  parsed  among  them 
equally. 
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D.  MATLAB  SPECIFICS 


1.  Sparse  Matrices 

The  lumped  mass  and  global  stiffness  matrices  are  sparse,  the  former  is  diagonal 
and  the  latter  is  block  tri-diagonal.  MATLAB  does  allow  assembly  of  a  full  stiffness  matrix 
followed  by  K= sparse  (K),  but  the  calculation  and  assembly  designed  from  the  outset  to 
be  sparse  is  a  more  efficient  use  of  computing  resources.  To  do  so,  the  global  indices  of 
each  element  of  the  collection  of  elemental  K  matrices  are  stored  in  vectors  i  v  and  j  v 
with  the  corresponding  values  stored  in  k  v.  After  the  elemental  calculations  are  complete, 
K=sparse  (i_v,  j_v,  k_v)  returns  the  sparse  global  stiffness  matrix. 

2.  CA  Implementation 

The  CA  portion  of  the  fluid  model  is  updated  using  sets  of  node  indices.  The  domain 
is  comprised  of  N  nodes  where  N  —  Nx  x  Ny  x  Nz.  In  this  work  Nx  —  Ny  —  Nz  =an 
odd  number;  this  keeps  the  eight  corners  of  the  domain  in  the  odd  set  and  is  convenient, 
not  required.  The  CA  coordinate  array  is  used  to  match  node  numbers  of  corresponding 
geometric  points  between  the  two  fluid  domains;  an  N  x  6  array  named  neighbor  tracks 
the  eponymous  relations  by  node  number  with  a  0  entry  indicating  the  end  of  the  domain 
in  that  direction,  the  velocity  potential  in  CA  domain  is  described  by  the  two  A-vectors, 
phi  and  phi  old.  A  large  portion  of  the  CA  set-up  is  the  definition  of  index  sets  for  these 
two  vectors.  The  largest  two  are  the  odd  and  even  interior  points,  oint  and  eint.  The  six 
faces  and  twelve  edges  of  the  rectangular  fluid  domain  are  likewise  split  into  odd  and  even 
index  sets.  This  infrastructure  makes  an  iteration  of  the  CA  rule  a  simple  matter  of  calling 
subroutines  UpdateFace  or  UpdateEdge  with  arguments  of  phiold,  neighbor, 
the  index  set  of  the  area  to  be  updated,  and  a  flag  indicating  the  boundary  condition  to  be 
employed. 
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E.  MATERIAL  PROPERTIES 

The  following  tables  specify  the  material  properties  used  for  various  calculations  in 
this  study. 


Property 

Value 

GXz 

345  MPa 

r< 

^Tyz 

345  MPa 

Table  2:  Material  Properties  of  aluminum  honeycomb.  From  [41],  [42] 


Property 

Value 

E 

72.4  GPa 

V 

0.3 

Table  3:  Material  Properties  of  aluminum  skins,  From  [41],  [42] 


Property 

Value 

Ex 

17.24  GPa 

P 

2020 

ITT3 

L'xy 

0.3 

Gxy 

6.619  GPa 

Ey 

17.24  GPa 

Ez 

7.929  GPa 

0.24 

Vyz 

0.24 

Gxz 

2.896  GPa 

r 

^Tyz 

2.896  GPa 

Table  4:  Material  Properties  of  E-glass,  From  [56] 
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Property 

Value 

E 

8.34  GPa 

P 

1180^1 

V 

0.28 

Table  5:  Material  Properties  of  Epoxy  Resin,  From  [57] 
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